{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimizing geometric parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pybamm.InputParameter is a powerful feature that enables efficient model reuse. It allows you to run multiple simulations with different parameter values without rebuilding the model, which significantly reduces computational overhead. However, it is not compatible with the standard meshes and finite volume code available in pybamm. This notebook shows how to use an alternative mesh implemented to allow for the use of pybamm.InputParameter. Unfortunately, this mesh is not compatible with the standard pybamm models, so custom models must be implemented for use. This limitation may be removed in future versions of pybamm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import pybamm\n",
    "\n",
    "domain_left = \"left domain\"\n",
    "domain_right = \"right domain\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model\n",
    "\n",
    "The model we will use is a one dimensional model of diffusion in cartesian coordinates. The model has two domains (left and right) with different diffusion coefficients, $D_l$ and $D_r$, and lengths. It is parameterized by the x coordinate of the left, right, and middle points, denoted $L_{min}$, $L_{mid}$, and $L_{max}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_left = pybamm.Variable(\"phi_left\", domain=domain_left)\n",
    "x_left = pybamm.SpatialVariable(\"x\", domain=domain_left, coord_sys=\"cartesian\")\n",
    "\n",
    "phi_right = pybamm.Variable(\"phi_right\", domain=domain_right)\n",
    "x_right = pybamm.SpatialVariable(\"x\", domain=domain_right, coord_sys=\"cartesian\")\n",
    "\n",
    "\n",
    "class TwoTankModel(pybamm.BaseModel):\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "        x = pybamm.concatenation(x_left, x_right)\n",
    "        D_left = pybamm.Parameter(\"D_l\")\n",
    "        D_right = pybamm.Parameter(\"D_r\")\n",
    "\n",
    "        left_flux = -D_left * pybamm.grad(phi_left)\n",
    "        right_flux = -D_right * pybamm.grad(phi_right)\n",
    "        flux = pybamm.concatenation(left_flux, right_flux)\n",
    "        source_left = pybamm.PrimaryBroadcast(pybamm.Scalar(10), domain_left)\n",
    "        source_right = pybamm.PrimaryBroadcast(pybamm.Scalar(0), domain_right)\n",
    "        source = pybamm.concatenation(source_left, source_right)\n",
    "\n",
    "        self.rhs = {\n",
    "            phi_left: -pybamm.div(left_flux) + source_left,\n",
    "            phi_right: -pybamm.div(right_flux) + source_right,\n",
    "        }\n",
    "        self.initial_conditions = {phi_left: 0, phi_right: 0}\n",
    "        self.boundary_conditions = {\n",
    "            phi_left: {\n",
    "                \"left\": (pybamm.Scalar(0), \"Dirichlet\"),\n",
    "                \"right\": (\n",
    "                    (D_right / D_left)\n",
    "                    * pybamm.boundary_gradient(phi_right, side=\"left\"),\n",
    "                    \"Neumann\",\n",
    "                ),\n",
    "            },\n",
    "            phi_right: {\n",
    "                \"left\": (pybamm.boundary_value(phi_left, side=\"right\"), \"Dirichlet\"),\n",
    "                \"right\": (pybamm.Scalar(0), \"Dirichlet\"),\n",
    "            },\n",
    "        }\n",
    "        phi = pybamm.concatenation(phi_left, phi_right)\n",
    "        my_grad_phi = pybamm.grad(phi)\n",
    "        self.variables = {\n",
    "            \"my_grad_phi\": my_grad_phi,\n",
    "            \"phi_left\": phi_left,\n",
    "            \"phi_right\": phi_right,\n",
    "            \"phi\": phi,\n",
    "            \"source_left\": source_left,\n",
    "            \"source_right\": source_right,\n",
    "            \"source\": source,\n",
    "        }\n",
    "\n",
    "        # define geometry\n",
    "        grad_phi_left = pybamm.grad(phi_left)\n",
    "        grad_phi_right = pybamm.grad(phi_right)\n",
    "        self.variables.update(\n",
    "            {\n",
    "                \"x\": x,\n",
    "                \"grad_phi_left\": grad_phi_left,\n",
    "                \"grad_phi_right\": grad_phi_right,\n",
    "                \"left_flux\": left_flux,\n",
    "                \"right_flux\": right_flux,\n",
    "                \"flux\": flux,\n",
    "            }\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-symbolic mesh\n",
    "\n",
    "We begin by showing how this process would work without using the symbolic mesh. We will use the default mesh, which is a uniform mesh with 50 points. We put the code in a function to show difference in time needed to change the parameters later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x346be9350>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGwCAYAAACOzu5xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAASJpJREFUeJzt3Qd4lFXCxfFDOoEkECChhCq9dwiIoiCKrmsDURFBBIW1uzbcz7oq9lXURcRVbEixN+wUlYSO9N57T0JCer7n3kkiKCqQTN53Zv6/55llJhsnlxgzZ+697z3lCgoKCgQAAOADgpweAAAAwIkiuAAAAJ9BcAEAAD6D4AIAAHwGwQUAAPgMggsAAPAZBBcAAOAzQuRH8vPztWPHDkVFRalcuXJODwcAAJwAc6RcWlqaatasqaCgoMAJLia01K5d2+lhAACAU7B161YlJCQETnAxMy1Ff/Ho6GinhwMAAE5AamqqnXgoeh0PmOBStDxkQgvBBQAA33Ii2zzYnAsAAHwGwQUAAPgMggsAAPAZfrXHBQDgveMmsrOznR4GfFRoaKiCg4NL5bkILgCAP2UCy8aNG214AU5VpUqVVL169RKfs0ZwAQD86cFgO3futO+WzeWqf3U4GHC8n6GMjAzt2bPHPq5Ro4ZKguACAPhDubm59kXHnGgaGRnp9HDgo8qXL2//NOElLi6uRMtGRGcAwB/Ky8uzf4aFhTk9FPi4yMLgm5OTU6LnIbgAAP4S/W9wy88QwQUAAPgMggsAAPAZBBcAAP5EvXr19Pzzzxc/3rVrl8455xxVqFDBXuIbyEs/H3/8cZl/XYILAp65VC8lI0eZOZ5NiAB8X8+ePXXbbbf97uMTJkw46bAxb948XX/99cWP//Of/9hLxBcvXqw1a9aUynhx4rgcGgEhKzdP2w8e0ZYDGdp6IKPwz18fp2XlKjwkSF0aVFHPxtXUs0k11a9agQ2JAFStWrVjHq9fv14dOnRQo0aNSnSoH1dqnRpmXOC3cvPy9d7cLTrrmRlqev9XOvvZmRryxjzd/8lyjf9xo75avksrdqba0GJk5eZr1pq9euTzFfZzz3x6hu7/eJm+X7lbGdmezwECnT1MLDvXkZv52qVtyJAhuvjii/XMM8/Yg9GqVKmiG2+88ZhLdo9eKjL3P/jgA7311lv2jY35540tW7booosuUsWKFRUdHa3LL79cu3fvLn6Ohx56SG3bttVrr72m+vXrKyIiwn7cPMe4ceP0t7/9zV4u3KxZMyUlJWndunV21sgsR3Xr1s2GpT+yadMm+zxTpkxRjx497JkpnTp1srNBZraoY8eOdlx9+/bV3r17j/lnzXjM1zTjadq0qf773/8eE65uuukm+30x/3/dunU1evToY/75ffv26ZJLLrFjN0Hu008/lbcx4wK/Y365/bBqj56Ytkpr9xwu/nhkWLBqV45U7dhI1bG38qpTJdJ+LKFypLYdzNCM1Xs1Y80ezd14wM7GvJ282d7CzGxM/Vj1aVFdV3aqrZBgMj8C05GcPDV/4GtHvvaKR85VZFjpv2xNnz7dvjibP01gGDBggA0Zw4cP/93nmiBwzTXX2HDywgsv2JBgqhCKQsvMmTPtoX0m/JjnmTFjRvE/a57bhJ4PP/zwmAPY/v3vf+u5556zt3vuuUdXXXWVGjRooFGjRqlOnToaOnSoDRDTpk3707/Hgw8+aANW0T9jnicqKsqO0wQLE6YeeOABjR071n7+u+++ax+/9NJLateunRYtWmT/ziYsDR48WGPGjLFBxAQi85xbt261t6M9/PDDeuqpp/T000/rxRdf1MCBA7V582bFxsbKWwgu8CtLth3S41+uVPKGA/ZxpchQ3Xx2I/29TU1VrRj2p0s/jeKj7G34GQ2UnpWr2ev3a8bqPTbMbD90RD+u3WdvHyzYpucHtFW9qhXK8G8GwFsqV65sX7xNmDCzDhdccIG+//774wYXs2wUHh5uA4vp3TG+/fZbLV261PY5mVoEw8zItGjRwgYdM/tRNINhPv7bpadrr73WhgrDBJfExETdf//9Ovfcc+3Hbr31Vvs5f+XOO+885p+58sor7d+je/fu9mPXXXed3eNzdNB59tlndemll9rHZiZoxYoVdgbIBBczi2RmUU4//XT7u9PMuPyWmXEyX8d4/PHHbdiZO3euzjvvPHkLwQV+wexTefrr1fr0lx32sZkhubZ7Pf2jZ0PFlA896eerEB6ic5rH25uZwVm/97C+W7lHL09fp8VbD+n8MT/qgb8114BOtdkHg4BSPjTYznw49bW9wQSMo2dAzOyLCSInauXKlTawFIUWo3nz5nYTsPn/ioKLeeH/bWgxWrduXXw/Pj7e/tmqVatjPpaZmanU1FQ70/NHTuR5ivqC0tPT7fKTCTNHBzQzWxQTE1McSszVU02aNLFBxCxn9enT5w+/ppmpMeMr+hreQnCBTzuUka2Xflint5I2KzsvXyZDXNK2lv55bhPVquTpxigpE0waxkXZ24VtauqfUxbbGZ17P1xqw8yTl7VSlYrhpfK1ALcz/z14Y7mmtJkX0JSUlN99/NChQ8UvzEVCQ0N/93f0RhO2eWE/nqO/ftEboeN97K/GFHoCz1P0HIcPe5bRx48fry5duhzzPEUhrn379nYWySxRfffdd3ZWqHfv3nr//feP+zV/+zW8xf0/fcAfmLZ0p+75YIlSMz0bZ7s3rKJRfZupZa1jfymVJhOGJg7rqtd+2mBneL5buVvnPn9IT/drrbOaxnnt6wI4OWaW4JtvvvndxxcuXKjGjRuX6tcym1uL9n8UzbqYJRcTkszMixvFx8fb4swNGzbYfSl/FgDNXh1z69evn515OXDggFf3sPwVggt80ocLt+nOqb8ov0BqWj1K9/ZtqjMbVyuTZZugoHK6/ozTdHrDarpt8iKt2X1Y106Yp0Fd6+q+85upfJh3prMBnLiRI0fafSu33HKLhg0bZvelfPHFF3rvvff02WeflerXMrMQZknGBACzOdYst/zjH//QmWeeaa/ocauHH37Yfn/MDJQJJFlZWZo/f74OHjyoO+64w24WNstmZuNuUFCQpk6davf1OH3oHpdGwOeYS5z/WRhaLu+YoC9u6aGeTeLKfK9J85rR+vSm0+1eGsNcfXTBiz/aDcIAnGWuypk1a5ZWrVplg4VZDjFXx5gX39LeOGp+93zyySd2k+8ZZ5xhv575+pMnT5abDRs2zF4O/cYbb9jgZYKW2bxrNuka5ookc8WQCV9mn4657PrLL7+0IcZJ5Qq8cWG8Q8zGJZMczbrmn21ggu96c/YmPfjpcnvfzHA8/PcWdgbEaT+u3WtngHanZikkqJyeG9DWXskE+DqzKdTsczj67BGgtH+WTub1mxkX+IxxM9cXh5bhPerrkYvcEVqMHo2q6evbztB5LaorN7/AbuD9ed0+p4cFAH6H4ALXM5OCY75fq9HTVtnHN5/d0O4lcdtlyJUiw/Tfge11Qesayskr0A1vL9DyHb+/qgEAcOoILnB9aDFX7zz3rafI7M4+jfXPPk1cF1qKmBmg5y5vo64NYnU4K9dWDJgzZgAApYPgAleHln9/vlL/neHp6Pi/C5rpprNPvdSsrISHBGvcoI72aqe9aVka/MZcHUzPdnpYQIn40XZI+PjPEMEFrpSfX6D/+3iZXv95o33874taaFiPBvIV5rTeCdd2Vs2YCG3Ym67r3pynI9l5Tg8LOGlFh5GZ4+qBksjIyDjuoXUni3Nc4EoPfLpM787ZYk/CffLS1rq8069HafuK6jERenNoZ/V7JUkLtxzSze8t0itXt6egET4lJCTEFvSZVmHzguP0pbDw0UbxjAxbBWDOgDm6XuFUcDk0XMf0Dd3y3iIbWv5zeVtd3K6WfNm8TQd09WtzlJWbrys719Hjl7R07R4d4HjMbIu5jNXbR7nDv1WqVMkeYHe8338n8/rNjAtcxWxk/deHnnKzm89q6POhxehUL1YvXNFOI99dYA/Pqx4doVt7u3+vDlAkLCzMtgSzXIRTZWbrSjrTUoTgAtfIzcvXrZMWKS0rVx3qVtYtvfznxf28ltX1yN9b6P5Plus/361RfHS4ruhcx+lhASfMLBFxAB3cgMVKuMYL36+1e0GiIkL0/IC2frcXZFBiPd141mn2/r8+Xqbpq7xb/Q4A/si/Xhngs5LW79dL09fZ+49f0kq1YyPlj+7s00T9OiQoz5yuO/UX7T+c5fSQAMCnEFzgOHPGye2TF8tsE+/fIUEX+nHHj9mUZoKZOePlQHq2Hvl8hdNDAgCfQnCBo8xFbfd8sES7UjPVoGoFPfT3FvJ3YSFBevKy1jI1S58s3qHvV+52ekgA4DMILnCUOavlmxW7FRpcTmOubKcK4YGxX7xN7UrFB+r966NlSs3McXpIAOATCC5wzJrdafp34VLJ3ec2VctaMQokt/durHpVIu1s0+gvPQWSAIA/R3CBIzJz8uwhc+ZQth6Nquq60+sr0JQPC9YTl7W29835LrPX73N6SADgegQXOGL0lyu1aleaqlYM07OXt7GtyoGoa4MqGtjFc57LvR8sVUZ2rtNDAgBXI7igzH23YrfeTNps7z/Tv43iogL7UKt7+za1ZYxbDmTouW/WOD0cAHA1ggvK1KGMbN39wRJ73ywP9WwSp0AXFRGqxy5pZe+bNuxFWw46PSQAcC2CC8rUmO/X2fNLGsdX1N3nNXF6OK5xVtM4XdKulvILpLvfX6Ks3DynhwQArkRwQZnZtC9dbydvsvf/74LmCg8pncItf/HA35rbPT9r9xzWy9PXOz0cAHAl1waXJ554wp4yettttzk9FJSSp75epZy8Ap3RuJq94ViVK4Tp4b+3tPf/O32dVu5MdXpIAOA6rgwu8+bN07hx49S6tedSUfi+BZsP6Mulu+xpsfed39Tp4bjW+a2qq0/zeOXmF9glI9OYDQBwcXA5fPiwBg4cqPHjx6ty5cpODweldKz/o1+stPcv71hbTatHOz0k1zKzjI9e3FLRESFauj1F//tpo9NDAgBXcV1wufHGG3XBBReod+/ef/m5WVlZSk1NPeYG9zEzLYu2HFL50GDdcU5jp4fjenHREfq/vzW39//z3RrtSct0ekgA4BquCi6TJk3SwoULNXr06BP6fPN5MTExxbfatWt7fYw4OebqmCe+8sy23HBmA/uijL9mWrLb1amkzJx8/ZeNugDgvuCydetW3XrrrXr33XcVEXFiL26jRo1SSkpK8c08B9zl7aTN2nrgiOKiwnX9GZ5SQZzYklHR7NTEOVu0M+WI00MCAFdwTXBZsGCB9uzZo/bt2yskJMTeZs6cqTFjxtj7eXm/P9ciPDxc0dHRx9zgrsPmxny/1t7/Z5/GigwLjObn0nJ6w6rqXC9W2Xn5eumHdU4PBwBcwTXBpVevXlq6dKkWL15cfOvYsaPdqGvuBwdz5oevefGHdUrNzFXT6lHq14FlvFOadenjmXWZMn+rth7IcHpIAOA417wFjoqKUsuWnjMsilSoUEFVqlT53cfhfpv3p+utJM9hc/ed30zBAVqiWBoljN0bVtHP6/brxR/W6ql+bZweEgA4yjUzLvAvT321msPmSskd53iqET5YuN2ePgwAgczVwWXGjBl6/vnnnR4GTuGwuS+W7uSwuVLSoW5lndWkmvLyC/RC4Z4hAAhUrg4u8D0cNufdWZePF2/Xuj1pTg8HABxDcEGp4rA572iVEGOrAAoKzKF0zLoACFwEF5TqYXNPfrXK3uewudJ3e2EQ/GLJTgoYAQQsggtKzfsLtmnLgQwOm/OSZjWidUHrGvb+f75d4/RwAMARBBeUivz8Ar32o6cQcMSZp3HYnJfc3ruR3fT8zYrdWrotxenhAECZI7igVHy3crc27ku3rcYDOnHYnLc0jIvSRW1r2fvPfbva6eEAQJkjuKBUFM22DOxaVxXCmW3xplt7NbIH+k1fvVcLNh90ejgAUKYILiixxVsPae6mAwoNLqch3eo5PRy/V69qBV3WnlkXAIGJ4IISG//jBvvn39vUUjxXEpWJm89uZIOiqQJI3rDf6eEAQJkhuKBETPHftKU77f1hPeo7PZyAUTs20h7wZzz3zRp78B8ABAKCC0rk9Z83Kr9A6tGoqr1cF2XnprMbKiwkyC7Tzd14wOnhAECZILjglKVk5GjyvK32/vAenNtS1mrElC/e6/JW0manhwMAZYLgglM2ce4WZWTnqWn1KDvjgrJ3TaJnM/RXy3dpV0qm08MBAK8juOCUZOfma8JszyXQw3o0ULly5ZweUkAyy3Od68fa5uiJc5h1AeD/CC44JZ/9skO7U7Ps8f5/b1PT6eEEtMGFsy5mBsz0RQGAPyO44KSZK1iKLoEe0r2e3SAK5/RpEa/46HDtO5ytr5btcno4AOBVvOLgpP20bp9W7UpTZFiwBnau6/RwAl5ocJAGdvH8e3hz9ianhwMAXkVwwUkbX3i8vzlHJCYy1OnhQNIVnWvbA+kWbjlE+SIAv0ZwwUlZtStVs9bstQ3F153OgXNuERcVofNb1bD330pi1gWA/yK44JTKFPu2rGFPb4X7Lo3+5JcdOpie7fRwAMArCC44YbtTM/XJ4u32Psf7u0/7OpXUsla0vVR98nzPwYAA4G8ILjhhZuNnTl6BOtWrrHZ1Kjs9HPyGOUunaNbl7aTN9mwXAPA3BBeckPSsXL07Z0vxgXNwJ3OmTqXIUG0/dEQ/rNrj9HAAoNQRXHBC3l+wTSlHclSvSqR6N4t3ejj4AxGhwRrQydMazSZdAP6I4IITOnDuvbme2ZZru9dXsLmkCK51dZe6Mg0MP67dp3V7Djs9HAAoVQQX/KVl21PtgXPmhNyL23raiOFe5mqvXk09s2LvJNNfBMC/EFzwl6YUXqFybovqHDjnIwZ3q1u8xHc4K9fp4QBAqSG44E9l5uQVXwLdv0OC08PBCep+WlU1qFbBhpaPFm5zejgAUGoILvhT367YrdTMXNWMiVD3hlWdHg5OUFBQOV3TtbC/KGmz3acEAP6A4II/NXWB5936ZR0S2JTrY8y/swphwXaDbtL6/U4PBwBKBcEFf2hnyhH9uHavvd+PZSKfExURqkvbe/69vcml0QD8BMEFf+jDhdtlVhg6149V3SoVnB4OTsE1iXWLl/zMoXQA4OsILjgusydiauHVRGzK9V2N4qPU7bQqMqf/P/fNGqeHAwAlRnDBcc3bdFCb9mfYPRLnt6rh9HBQAv/s08QeSPfBwm2aTg0AAB9HcMFxFc22XNC6hiqEhzg9HJRAh7qVdV13T5v3qA+X2uoGAPBVBBcct1Dxi6U77f3+HT29N/D9WZf6VStoV2qmHvtihdPDAYBTRnDB73y5dKcysvPsC13HupWdHg5KQfmwYD3Vr7VdMpoyf5tmrGbJCIBvIrjgd6bO31Z8CXQ580oHv9CpXqyu7fbrklFqJktGAHwPwQXH2LQvXXM3HZA5a+7S9hQq+pu7zm2iulUitTMlU49/sdLp4QDASSO44BimlM/o0aiaasSUd3o48MKS0dP92tglo0nztmrWGs8BgwDgKwguKJaXX2AvmTX6d+TsFn9lDhQcnFjP3r/3gyVKY8kIgA8huKDYz+v22SWEmPKh6t0s3unhwIvuPq+J6sRGaodZMvpyldPDAYATRnDB7woVL2pbUxGhwU4PB14UGRZirzIy3pu7RT+t3ef0kADghBBcYKVk5Ojr5bvs/f4dOLslEHRtUEWDC7uM7vlgiQ5n5To9JAD4SwQXWJ/+sl3ZuflqWj1KLWtFOz0clJG7z2uq2rHlbQHj6C+5ygiA+xFccMwykTkpl7NbAoepc3jyMs+S0btztth9TgDgZgQXaPWuNC3ZlqKQoHK6uG1Np4eDMtbttKoa1NWzZHT3+ywZAXA3gguKCxV7NYtTlYrhTg8HDri3b1MlVPYsGT05jauMALgXwSXAmbNbPvllh73PptzAXjJ6qnDJ6O3kzZq9niUjAO5EcAlwCzYf1N60LEVHhOiMxtWcHg4c1K1hVQ3sUsfev2vqEm09kOH0kADgdwguAc40QRvnNK+usBB+HALdqPObqV6VSLtk1O+V2VqzO83pIQHAMXilCmD5+QX6apnn7JbzW1V3ejhwgYrhIZp8Q6Iax1fU7tQsXT4uSYu3HnJ6WABQjOASwBZtPaRdqZn2xer0RlWdHg5cIj46QlNuSFTb2pV0KCNHV41P5jJpAK5BcAlg0wqXiczVROEhHPGPX1WKDNO7w7ro9IZVlZGdp2vfmFc8OwcATiK4BKiCggJNK3wh6tuyhtPDgUuvNPrfkI7q27K6svPy9Y93F2jKPM+l8wDgFIJLgDIHzpkNmJFhwerZhKuJcHxmJu7FK9tpQMfayi+Q7v5giV77cYPTwwIQwAguAerLZZ5lorOaxtEEjT8VEhykJy5rpevPaGAfP/rFSj3z9Wo7awcAZY3gEoDMC07x1UQsE+EEmP6qUX2b6u7zmtjHL01fp/s/WWavTAOAskRwCUArdqZq8/4MRYQGsUyEkwov/+jZUI9d0lKmh/Od5C26dfJi2yoOAGWF4BKApi31zLb0bBxnN2ACJ2Ngl7oac0U7W8r52S87dP3b83UkO8/pYQEIEASXAFwmKjotty+HzuEUXdimpl4b3NHO2s1YvVeD/jdHKUdynB4WgABAcAkwa3Yf1oZ96QoLDtLZTeOcHg58WM8mcXrnui6KigjR/M0HdcWrybb3CgC8ieASYIpmW85oXFVREaFODwc+rmO9WE2+PlFVK4Zr5c5U9X9lNuWMALyK4BJgphVeBs2hcygtzWtG6/0RiUqoXF6b9meo/ytJWks5IwAvIbgEkHV70uxSUWhwOfVuFu/0cOBH6lWtoPdHdFOjuIq2/6o/5YwAvITgEoBXE3VvWFUxkSwToXRVj/GUM7YpLGccOD5ZsylnBODPwWXs2LFq3bq1oqOj7S0xMVHTpk1zelh+o6ibiEPn4C2VK3jKGbs3rKL07DwNeWOevl5OOSMAPw0uCQkJeuKJJ7RgwQLNnz9fZ599ti666CItX77c6aH5vE370u3Bc8FB5XROc5aJ4D0Vw0P0+pBOOrdFvC1nHPnOAk2dTzkjAD8MLhdeeKHOP/98NWrUSI0bN9Zjjz2mihUrKjk52emh+c1sS7fTqth3xYC3yxlfvqq9+ndIsOWMd71POSOA0uHaY1Pz8vI0depUpaen2yWj48nKyrK3IqmpqWU4Qt/C1URwopzxqX6tFVM+VK/9tNGWM5pD6u44p7GtDwAAn59xMZYuXWpnWcLDwzVixAh99NFHat68+XE/d/To0YqJiSm+1a5du8zH6wvMuRpLtqUoqJzUpwXLRCg7JqD864JmuutcTznjiz+s0wOfLKecEYD/BJcmTZpo8eLFmjNnjkaOHKnBgwdrxYoVx/3cUaNGKSUlpfi2dSvr6MdT1ATduX6sPSgMKOvwcuNZDfXviz3ljG8nb9btUxYrJ49yRgB+sFQUFhamhg0b2vsdOnTQvHnz9MILL2jcuHG/+1wzK2Nu+HNfFi4Tnd+KZSI4Z1DXunbZ6I7Ji/XJ4h1Ky8y1+2DKhwU7PTQAPsR1My6/lZ+ff8w+FpycHYeOaNGWQ/ad7rktKFWEs/7epqbGX+MpZ/xh1R4Nfn2uUjMpZwTgo8HFLP3MmjVLmzZtsntdzOMZM2Zo4MCBTg/N55eJOtatrPjoCKeHA+ispnF6u7Ccce6mA7piHOWMAHw0uOzZs0fXXHON3efSq1cvu0z09ddf65xzznF6aD6Lq4ngRp3qxWrS9V1VtWKYPV/o8nFJ2naQckYAf61cQUGB32zvN5dDm6uLzEZdc/JuoDPvYjs//p3Mv+HZ956tmpXKOz0k4Bgb96Xr6tfmaPuhI6oeHaF3hnVWw7gop4cFwMWv366acUHpmrlmrw0tLWtFE1rgSvVNOePIRDUsKmd8JUlLtlHOCOCPEVz82PTVe+yfZzeJc3oowB+qEVPeU86YEKODGTm68tVkzV5POSOA4yO4+KncvHzNWrPX3j+T4AKXizXljMO72kqKonLGbyhnBHAcBBc/tXDLIXtORqXIULWtXcnp4QAnXM7Yp3m8snPzNfLdhfpgwTanhwXAZQgufmpG4TLRmY2r2UZowBdEhAbrvwPbq1+HBOXlF+ifU3/R6z9tdHpYAFyE4OKnpq/2LBOdxTIRfLGc8bLWuu70+vbxI5+v0HPfrpEfXQAJoAQILn5oV0qmVu5MtaflntG4mtPDAU5aUFA5/d8FzXRnn8b28Zjv1+qhTylnBEBw8etlojYJleymR8BXyxlvOruR/n1RCxvC30zarDsoZwQCHsHFjy+DZpkI/mBQYj09P6CtQoLK6ePFO3TD2wuUmZPn9LAAOITg4mfM1Rg/r9tv75/VlGUi+IeL2tbSq9d0UHiIp5zxGsoZgYBFcPEz8zcf0OGsXNsB07JmjNPDAUrN2U3jPeWM4SGau/GAPahu32HKGYFAQ3DxMzMKryY6s3Gc3eAI+JPO9WP13vVdVaVCmJbvSNXlryTZniMAgYPg4memryrc38IyEfxUy1oxmjoiUbUqldeGfenqN3a21u057PSwAJQRgosf2XYwQ2v3HJaZaOnRkOAC/9WgWkUbXk6rVkE7UzJ1+bgkLd2W4vSwAJQBgosfLhN1qFtZMZGhTg8H8CrTeD51RDe1qhWjA+nZunJ8spLWezamA/BfBBc/PL+lJ5dBI0CYc4omDu+irg1i7ab0wW/M1bcrdjs9LABeRHDxE+Zci+LLoAkuCCBREaGacG1n9W7mKWcc8c4CfbiQckbAXxFc/IS5PPRITp7io8PVrEaU08MByryc8ZWr2+vS9rVsOeMdU37RhJ8pZwT8EcHFz/a3mNkWc1Q6EIjljM/0a6Mh3erZxw99tkLPf0c5I+BvCC5+t7+Fq4kQuMzZRQ9e2Fy39/aUMz7/3Vo9/NkKyhkBP0Jw8QOb9qXb8yxMl0v3hlWdHg7gKDPjeGvvRnrowub28YTZm3Tn1F8oZwT8BMHFj2ZbOtWLtRsVAUhDutfXfwa0UXBQOX24aLtGvrOQckbADxBc/MD0ov0tnJYLHOOSdgkad7WnnPG7lbs15I25SqOcEfBpBBcfdyQ7T8kbuAwa+CO9m8frzaGdVTE8RMkbDuiq8XO0n3JGwGcRXHycCS1Zufm2t6VhXEWnhwO4UtcGVTTp+q72wLql21NsRcAOyhkBn0Rw8XHTj7qaiMuggT8vZ5xyQ6JqxkRo/V5POeOGvZQzAr6G4OLDzPkUPxS1QbNMBPwlMys5dWQ3NahaQTtSMtX/lSQt2045I+BLCC4+zLxr3HbwiMKCg9StYRWnhwP4BLOsOmVEolrWitZ+U874arLmFO4TA+B+BBc/uAy6S4NYRYaFOD0cwGdUrRiu94Z3Vef6sUrLytU1r8/VD6soZwR8AcHFT475B3ByzJlHbw3trF5N4+wG9+vfWqCPF213elgA/gLBxUelZ+VqzkbP9DbH/AMlKGcc1EGXtKul3PwC3TZ5sd5K2uT0sAD8CYKLj5q9fr9y8gpUt0qk6let4PRwAJ8VGhykZ/v/Ws74wCfLNeb7tZQzAi5FcPFRP6/bZ//s0agql0EDpVTOeGuvRvbxc9+u0b8/X0k5I+BCBBcflbTes0zU7TRKFYHSYN4A3H5OYz3wN0854+s/b9Rd7y9RLuWMgKsQXHzQvsNZWr07rfhEUAClZ+jp9e3SkSln/GDhNo18l3JGwE0ILj6oqJuoWY1oe4Q5gNJ1WYcEvXJ1B4WFBOnbFbt17RvzdDgr1+lhASC4+O7GXKPbacy2AN5yjilnvNZTzpi0Yb+uGp+sA+nZTg8LCHgEF5/e30JwAbwp8bQqmji8iypHhmrJNk85484UyhkBJxFcfIxptN24L92uv5tTPwF4V+uESpo6optqxERo3Z7D6jc2iXJGwEEEFx+dbWlVK8ae/AmgjMoZRyTaM5O2HzpiZ16W76CcEXACwcXHsL8FcEZC5UgbXprXiNa+w9m6Ylyy5m484PSwgIBDcPEh5iTPpPWeg+c4vwVwppxx0g1d1bmep5xx0P/maPoqT9kpgLJBcPEhm/dnaEdKpkKDy6lD3cpODwcISNERoXpzaGedXVjOOPyt+fpkMeWMQFkhuPjgMlG7OpVVPizY6eEAAcv89zduUAdd1LZmcTnj28mbnR4WEBAILj7EnCVhsL8FcEc5438ub6trEuvK9DHe//EyvfQD5YyAtxFcfAT7WwB3ljM+/PcWuuXshvbxM9+s0WNfrCS8AF5EcPERa/cctlcyRIQGqW3tSk4PB8BR5Yx39Gmi+wvLGV/7aaPuppwR8JqQkvzDS5Ys0Y8//qiwsDB169ZNLVq0KL2R4Riz13lmWzrVi7X9KQDc5brT6yumfKju+WCJpi7YptTMHL1wRTtFhLIfDXBFcHnhhRd0++23Kzo6WsHBwTp48KBatWqlN998U23bti3VQeLo81tYJgLcql+HBEVFhOjmiYv09fLdGjphnl69pqPtOwJQOk7qrfvrr7+uhQsXKisrS4899pieeOIJG1j279+vDRs2qG/fvurRo4dmz55dSsODkZdfUNwIzcZcwN3ObVFdE67tpAphwfYNx8DxyTpIOSPgTHB55pln1KVLF1WsWNGGlXnz5tmZl5kzZ6py5co2yJjbnXfeWXojhFbsSFVqZq6iwkPUoma008MB8Be6NayqicO72nLGXwrLGXelZDo9LCDwgsuKFSuUlpZmZ1RCQ0MVFBSkSZMm6fzzz1dsbKwaNGigjz76SAsWLNAXX3yhTZs2eW/kAWR24dVEXRrEKiSY/S2AL2hT25QzJqp6dITdXH/Z2Nm2IBVAyZz0q2BERIQ6deqk7t27q02bNkpOTrZhZunSpXr00UfVsGFD5eTk6JprrrFBxuyBQensb0lkfwvgUxrGRen9kb+WM/Z/JcnOoAI4daf89v3ZZ5/VU089pWHDhtl9L40aNdKFF15ol5Fq1qxpl5K2bNmiKVOmlGB4yMnL17xNniI39rcAvlnOOOWGonLGLA14NUnzC/+bBlCGwcVcOWSWhDZv3qyuXbvamZhKlSrpxRdf1JNPPmk/JyEhQeedd96pfgmYS863HVJGdp5iK4SpSXyU08MBcAqqRYXrveu7qlO9ykrLzNXVppxxNeWMwKko0YaJ0047Td9++622b9+uDz/80O53Wb9+vQYOHFiSp8VRZq8rXCZqUMWe0gnAN5kzXt4a2kU9m1RTZk6+hr85X5/+ssPpYQE+p1QOF4iPj9dFF11UGk+FP9zfwjIR4A/ljK8O6qg7p/5iQ8utkxYp9UiOru5a1+mhAT6DS1RcLDMnTwu2HLT32d8C+Adz8vXzA9pqUFdPOeP/fbxML09fR78RcIIILi62cPNBZefm28spzVUJAPyDWfZ95KIWurmwnPHpr1fr8S8pZwROBMHFR5aJTJEbAP9h/pv+Z58m+r8LmtnH43/caHuOKGcE/hzBxQcOnmN/C+C/hvVooKf7tZbZez9l/jbdNHGRsnLznB4W4FoEF5c6nJVrjwo32N8C+Lf+HWtr7NUdFBYcpK+W77LljOlZuU4PC3AlgotLzdt4wJYr1omNtAdYAfD/csY3CssZf163X1e9NodyRuA4CC4uXyZitgUIHN2PLmfcesiesrs7lXJG4GgEF5dK2sD5LUCgljOaigBzNeGa3Z5yxk2UMwLFCC4udCgjW8sLi9gILkDgaRQfZZul61WJ1LaDR9TvlSSt3Ek5I+C64DJ69GjbPB0VFaW4uDhdfPHFWr16tQJN8oYD9mCqRnEVFRcV4fRwADigdmykpo7opmZF5YzjkrRgM+WMgKuCy8yZM3XjjTcqOTnZdiDl5OSoT58+Sk8PrGnSZJaJABSWM066vqs61q2s1MxcDXxtjmau2ev0sABHuSq4fPXVVxoyZIhatGihNm3aaMKECdqyZYttoQ4kczd63lV1qU9wAQKdKWd8+7pfyxmHvTlPny+hnBGBy1XB5bdSUjznmMTGxh73/8/KylJqauoxN1+XlpmjVbs8f4+O9So7PRwALipn/FvrGsrJK9DN7y3SxDlbnB4W4AjXBpf8/Hzddttt6t69u1q2bPmHe2JiYmKKb7Vr15avW7TlkPILzPp2ecVHs78FwK/ljC9c0U4Du9Sxe+Du+2ipxs5Y7/SwgDLn2uBi9rosW7ZMkyZN+sPPGTVqlJ2VKbpt3bpVvm7+Zk8bdMe6x59lAhC4goPK6dGLW+rGs06zj5/8apVGT6OcEYElRC5000036fPPP9esWbOUkJDwh58XHh5ub/5k/ibP/haWiQD8UTnjXec2VaXyYXrsy5UaN3ODUjJy9NglrWywAfydq2ZczLsGE1o++ugj/fDDD6pfv74CSU5evhZvPWTvM+MC4M8MP6OBnrrMU844ad5W3fIe5YwIDCFuWx6aOHGiPvnkE3uWy65du+zHzf6V8uXLy9+ZA6YysvMUHRFiz3ABgD9zeafaii4folveW6wvlu5UamaOxg3qoMgwV/1qB/x3xmXs2LF2r0rPnj1Vo0aN4tvkyZMVCOZv8uxv6VC3soKY8gVwAs5rWUOvD+mkyLBg/bh2nz3rxZy+Dfgr1y0VHe9mznYJBPMLT8XsWI9lIgAn7vRGVfXusC72zBdzZeKAccnaQzkj/JSrgksgMwGtaMbFnJIJACejXZ3KtpwxLipcq3en2X6jLfsznB4WUOoILi6x9cAR7UnLUmhwOdsOCwAnq0n1KH0wspvqVonUlgMZ6vfK7OIDLQF/QXBx2TJRy1oxiggNdno4AHy5nPGGRDWtHmXfDJllo4VbPLO5gD8guLjEvMJlok7sbwFQQnHREZp8faLa16mklCM5Gjh+jn5cSzkj/APBxSWK6urNFUUAUFIxkaF6Z1gXndG4mo7k5GnohHn6culOp4cFlBjBxQXMpYtrdh+299mYC6C0mPNcXrumoy4oLGe8aeJCTZpLOSN8G8HFBYrWnxtUraAqFf2rwgCA8+WMY65opys717EFrvd+uFTjZlLOCN9FcHHR/hb6iQB4g+kwevySlhrZ01POOHraKlvQSDkjfBHBxQUWFJ/fwsZcAN4rZ7znvKYa1bepfTx2xnrd99Ey5ZlpGMCHEFwcZkrRFm/zFCt2YMYFgJfdcOZpeuLSVrac8b25W3TLpEXKzs13eljACSO4OGzZ9lT7SyO2Qpjd4wIA3nZF5zp66ar29sDLL5bs1LC35isjO9fpYQEnhODisPmbfr0M2kzlAkBZOL+Vp5yxfGiwZq3Zq6tfm6OUjBynhwX8JYKLw+ZvLjp4jmUiAGWrR6Nqene4p5xxoSlnfDWJcka4HsHFQWZH/4LC4NKBjbkAHNC+TmVNvqGrqkWFa9WuNPUfl6StByhnhHsRXBy0YV+6DqRnKzwkSC1rRTs9HAABqmn1aH0wopvqxEZq8/4MXTZ2tlbvSnN6WMBxEVxcsL+lTUIlhYdQrAjAOXWqROr9EYlqEu8pZ7x8XBLljHAlgouD5nPwHAC3lTPe0FXtCssZzYbdn9buc3pYwDEILi7YmEtwAeAWlSLD9O6wLurRqKoysj3ljNMoZ4SLEFwcsu9wljbuS7f3O9RhYy4Al5UzDu6oC1rVUHZevm6cuFBT5m11eliARXBxeJmocXxFWz8PAG5i9t2NubKdruhU25Yz3v3BEr06i3JGOI/g4pAFmz0bczvWY7YFgHvLGUdf2ko3nNnAPn78y1V6inJGOIzg4nQjdF32twBwL3Oi96i+zWxBo/HfGev1r48pZ4RzCC4OOJKdp+U7Uuz9Tsy4APABI3ueZmdfTDPJxDlbdCvljHAIwcUBv2w7pJy8AsVFhSuhcnmnhwMAJ+TKznX04pXtbDnj50t2ajjljHAAwcUBRcf8m9kWihUB+JK/ta6p1wZ7yhlnrtmrQf+bSzkjyhTBxQHzjmqEBgBfc2bjanpnWGdFR4TYN2K2nDGNckaUDYJLGcvP/7VYkf0tAHyVKYadfEOiqlYsLGd8hXJGlA2CSxlbsydNaZm5igwLVrMaUU4PBwBOWbMa0fpgZKJqx5a35Yz9XpmtNbspZ4R3EVwcOnjOdIGEBPPtB+Db6lapoPdHdLOHae5O9ZQzLqKcEV7EK6dDjdBmmhUA/EF8dISm3JCotrUr6VBGjgZSzggvIriUsQWF70Q4eA6AP5Yznt7w13LGr5ZRzojSR3ApQ/sPZ2nrgSP2fpvalZweDgCUqgrhIfrfkI7q27K6LWf8x7uUM6L0EVzK0JJtntNyT6tWQTHlKVYE4J/ljOaQugEdfy1nHD9rg9PDgh8huJShRVsP2T+ZbQHgz8yFB09c1krXn+EpZ3zsy5V6+mvKGVE6CC5l6JfC4GI2sAGA/5czNtXd5zWxj1+evl7/RzkjSgHBpYyYdxqmo8hok0BwARAY4eUfPRvqsUta2nLGd+ds0W2TF1POiBIhuJSRLQcy7GWCYcFB9tAmAAgUA7vU1Zgr2ikkqJw++2WHrn97vo5k5zk9LPgogksZWVy4TNS8ZrTCQvi2AwgsF7Yx5YwdFREapBmrTTnjHKUcoZwRJ49X0DIOLuxvARCoejaJ0zvXdVFURIjmbz6oK15N1t60LKeHBR9DcCnjjbltasc4PRQAcEzHerGafL2nnHHlzlT1f2U25Yw4KQSXMpCTl69lO1Lt/ba1OTEXQGAzS+bvj0hUQuXy2rQ/wzZLr6WcESeI4FIGVu1Ms7vooyNCVK9KpNPDAQDH1avqKWdsFFdRu1Iz1X9cUvGSOvBnCC5lYHHRZdC1K9nLAwEAUvUYTzljm6JyxvHJmr2Ockb8OYJLGeDgOQA4vsoVPOWM3RtWUXp2noa8MU9fL9/l9LDgYgSXMkBwAYA/VjE8RK8P6aRzW8TbcsaR7yzQ1PmUM+L4CC5elpaZo3V7D9v7rTkxFwD+sJzx5avaq3+HBFvOeNf7S/Taj5Qz4vcILl62dFuKTK9YrUrlVS0q3OnhAICryxmf6tdaw06vbx8/+sVKPfvNasoZcQyCSxltzGWZCAD+mrmA4V8XNNNd53rKGV/8YZ0e+GS58ilnRCGCi5dx8BwAnHx4ufGshvr3xZ5yxreTN+v2KYvtmVgAwaXMjvrn4DkAOBmDutbVC4XljJ8s3qHr36KcEQQXr9qVkqndqVkKKie1rEUjNACcrL+3qanx13jKGaev3qvBr89VaibljIGM4FIGsy2N46MUGRbi9HAAwCed1TRObxeWM87ddEBXjKOcMZARXLzoFzbmAkCp6FQvVpOu76qqFcO0YmeqLh+XpG0HKWcMRAQXL+LgOQAoPS1qxmjqiG72eImN+9JtOeO6PZQzBhqCi5fk5RdoybYUe9/0cAAASq6+KWccmaiGcRW1MyXThpclhbPbCAwEFy/ZsPewDmflqnxosG0/BQCUjhox5T3ljAkxOpiRoytfTdbs9ZQzBgqCi5c35raqFWNPgwQAlJ5YU844vKu6nfZrOeM3lDMGBF5Rvb0xtw7LRADgzXLGPs3jlZ2br5HvLtQHC7Y5PSx4GcHFyzMubShWBACviQgN1n8Httdl7RPs3sJ/Tv1Fr/+00elhwYsILl6QmZOnVTs9O9056h8AvMssxz/dr7WGdveUMz7y+Qo99+0ayhn9FMHFC5bvSFVufoE9b8BctgcA8K6goHK6/2/N9M9zGtvHY75fq4c/W0E5ox8iuHj5/BZTFgYA8D7z+/bmXo30yEUt7OMJszfZpSPKGf0LwcUL2N8CAM65JrGenh/QVsFB5fTRou0a8fYCu4QP/0Bw8eIVRRw8BwDOuLhdLb06qIPCQ4L0/ao9uoZyRr9BcCllB9OztXm/pz+DGRcAcE6vZvF6a2hnRYWHaO7GA7pqfLL2H6ac0dcRXLw029KgagXFRIY6PRwACGhdGlTRe9d3VZUKYVq2PVX9xyVp+6EjTg8LJUBw8db+FpaJAMAVWtYy5YyJqhkToQ1709V/7Gyt33vY6WHBH4LLrFmzdOGFF6pmzZp2d/jHH38sX72iyHRoAADcoUG1inp/ZDedVq2CdhSWMy7b7inChW9xVXBJT09XmzZt9PLLL8sXmcOOfqERGgBcqWYlTzmj6ZA7kJ6tK15NVvKG/U4PCycpRC7St29feztRWVlZ9lYkNTVVTtp28Ij9jyE0uJya14x2dCwAgN+rUjFcE4d30bA352vOxgMa/PpcvXxVe/VuHu/00OCLMy4na/To0YqJiSm+1a5d29HxLCpcJmpeI1rhIcGOjgUAcHxREaF6c2hn9W4Wr6zcfN3wzgJ9tIhyRl/h08Fl1KhRSklJKb5t3brVHftbWCYCANeXM75ydXtd2q6WLWe8ffIvmvAz5Yy+wFVLRScrPDzc3ty3MZfgAgC+UM74TP82ii4fausBHvpshQ4dydGtvRpR1+JiPj3j4iamC2PZDs/G3LZ1CC4A4CvljA9e2Fy39/aUMz7/HeWMbkdwKSXr9hxWZk6+PaGxfpUKTg8HAHCCzOzKrb0b6aELm9vHZvblTsoZXctVS0WHDx/WunXrih9v3LhRixcvVmxsrOrUqSM3KzoPwFxNZBI8AMC3DOle3554fufUJfpw0XalZubqpava2f0wcA9XzbjMnz9f7dq1szfjjjvusPcfeOABud3yHanFJzQCAHzTJe0SNO7qDgoLCdJ3K3dryBtzlUY5o6u4Krj07NnTHuL229uECRPkdkUzLi1rcX4LAPgyc6aLKWesGB6i5A2mnHEO5Ywu4qrg4qvMpXQrdhbOuNRkxgUAfF1XU844vKtiK4Rp6fYUXT4uSTsoZ3QFgksp2LQ/XRnZeYoIDbJ9GAAA39cqIcZWBJhyxvV709Vv7GxtoJzRcQSXUlwmalYjWsFszAUAv9EwrqKmjuymBlUpZ3QLgktpbsxlmQgA/E4tU844ItHuYdyfnq0rX03WHMoZHUNwKQVszAUA/1a1Yrjd89K5fqzSsnJ1zetz9cOq3U4PKyARXErIXPVUFFxaMOMCAH5dzmiuNurVNM6WM17/1gJ9vGi708MKOASXEtp28Ig9pCg0uJwax0c5PRwAgLfLGQd10CXtaik3v0C3TV6st5I2OT2sgEJwKaHlhf1ETapH2QOLAAD+LTQ4SM/2b6Mh3erZxw98slxjvl9rZ+DhfbzSltCy7Z6NuS1qsEwEAIFWzmiapI3nvl2jf3++knLGMkBwKaGiRmg25gJA4JUz3n5OYz3wN0854+s/b9Rd7y9RLuWMXkVwKa2NuXQUAUBAGnp6fbt0ZM7x+mDhNo18d6Eyc/KcHpbfIriUwJ60LO07nC1z5lyz6sy4AECguqxDgl4pLGf8dsVuXfvGPB3OynV6WH6J4FICRbMt5mTF8mHUngNAIDunebzevNZTzpi0Yb+uGp+sA+nZTg/L7xBcSoATcwEAR0s8rYomDu+iypGhWrLNU864M4VyxtJEcCmFGZfmNVkmAgB4tE6opKkjElUjJkLr9hxWv7FJlDOWIoJLacy4sDEXAHCUhnFRNrzUr1pB2w8dsTMvRed+oWQILqfIrFuaH0aDGRcAwG8lVI604aV5jWh7IccV45I1d+MBp4fl8wgup6goOderEqnoiFCnhwMAcGk546QbuqpzPU8546D/zdH0VXucHpZPI7iU9MRclokAAH/CvLl9c2hnnV1Yzjj8rfn6ZDHljKeK4FLCGZcWLBMBAP6COTJj3KAOuqhtzeJyxreTNzs9LJ9EcDlFXAoNADjZcsb/XN5W1yTWleljvP/jZXrpB8oZTxbB5RSkZeZo4750e58ZFwDAyZQzPvz3Frrl7Ib28TPfrNFjX6wkvJwEgsspWFE421IzJkJVKoY7PRwAgI+VM97Rp4nuLyxnfO2njbqbcsYTRnA5BcsKgwsbcwEAp+q60+vrmcJyxqkLtunGiZQzngiCyylYXnhiLvtbAAAl0a9Dgv47sL3CgoP09fLdGjqBcsa/QnApwcZc9rcAAErq3BbVNeHaTqoQFqzZ6/dr4PhkHaSc8Q8RXE7Skew8rd2TZu9z1D8AoDR0a1hVE4d3teWMvxSWM+5KyXR6WK5EcDlJq3alKr/AnIYYpvhoNuYCAEpHm9qVNOWGRFWPjtDaPYd12djZ2lR4BSt+RXA51Y25NWPsznAAAEpLo/govT8y0dbJmD68fq8kFV/JCg+Cy6luzK3F/hYAgLfKGbsVljNmacCrSZq/iXLGIgSXk7Ss8Kh/rigCAHhLtahwvXd9V3WqV1lpmbm6+n9zNGM15YwGweUkZOfma82uw8VLRQAAeEtM+VC9NbSLejappsycfA17c74++2WHAh3B5SSYq4my8/IVFRGi2rHlnR4OACAAyhlfHdRRF7bxlDPeMmmR3p0T2OWMBJeTsHz7r8WKbMwFAJSFsJAgPT+grQZ2qWPLGf/10TK9PH1dwPYbEVxOZX8LG3MBAGXI1AI8enFL3XSWp5zx6a9Xa/S0VQEZXgguJ2FZ8RVF7G8BAJQtM9N/57lN9K/zm9nHr87aoHs/WKo8c7hYACG4nCDzg7Fyp+fEXDbmAgCcMvyMBnrqstYKKidNnr9VN01cqKzcwClnJLicoI37DutITp7KhwarftUKTg8HABDALu9Uu7iccdqyXbpuwnylB0g5I8HlBC0r3JjbvGa0XWsEAMBJ57WsodeHdFJkWLB+WrdPA1+bo0MZ/l/OSHA52f0tNEIDAFzi9EZV9e6wLvbMl8VbD9lyxt2p/l3OSHA5ySuKWrAxFwDgIu3qVNbUEYm2+HfN7sPq98psbd7vv+WMBJcTkJ9fcMwZLgAAuEljU844opvqVonU1gOecsaVO/2znJHgcgK2HsxQWlau3QTVKL6i08MBAOB3aseacsZENa0epb1pWRowLkkLNvtfOSPB5QQsL6wUb1I9SqHBfMsAAO4UFxWhydcnqkPdyko15YyvzdXMNXvlT3gVPgGmm+isJtXsJigAANwsJjJUb1/XWWc2rmaP8Rj25jx9sWSn/EW5Aj86Lzg1NVUxMTFKSUlRdDRX/wAAAld2br7umLJYny/ZKVOv9/glrXRl5zry9ddvZlwAAPDTcsYXrminqwrLGUd9uFRjZ6yXryO4AADgp4KDyumxi1vqHz1Ps4+f/GqVRk9b6dPljAQXAAD8vJzx7vOaalTfpvbxuJkb7OyLr5YzElwAAAgAN5x5mp64tJUtZ5w0b6tufs83yxkJLgAABIgrOtfRy1e1V2hwOX25dJeGvTlfGdm+Vc5IcAEAIID0bfVrOeOPa/fpah8rZyS4AAAQYHo0qqZ3CssZF245pAHjkrXHR8oZCS4AAASg9nUqa8oNiYqLCtfq3Wm232jL/gy5HcEFAIAA1aS6p5yxTmykthzIsM3Sq3a5u5yR4AIAQACrUyVS7xeWM+6x5YzJWrjloNyK4AIAQICLi/aUM7avU0kpR3I0cPwc/bjWneWMBBcAACBTzmg27PZoVNWWMw6dME9fLnVfOSPBBQAAWJFhIXptcEdd0KqGcvIKdNPEhZo0d4vchOACAACKhYcEa8yV7XRl59oyrQD3frhU42a6p5yR4AIAAH5Xzvj4Ja004kxPOePoaatsQaMbyhkJLgAA4LjljPf2bWpvxtgZ63XfR8scL2ckuAAAgD9kZl1GX9pK5cpJ783dolsmLVK+g+GF4AIAAP7UlZ3r6KUrPeWMTeKjFGQqph0S4thXBgAAPuOC1jXsSbunVavg6DgILgAA4IQ0jKsop7lyqejll19WvXr1FBERoS5dumju3LlODwkAALiA64LL5MmTdccdd+jBBx/UwoUL1aZNG5177rnas2eP00MDAAAOK1fghouyj2JmWDp16qSXXnrJPs7Pz1ft2rV1880369577z3mc7OysuytSGpqqv3clJQURUdHl/nYAQDAyTOv3zExMSf0+u2qGZfs7GwtWLBAvXv3Lv5YUFCQfZyUlPS7zx89erT9ixbdTGgBAAD+y1XBZd++fcrLy1N8fPwxHzePd+3a9bvPHzVqlE1nRbetW7eW4WgBAEBZ8+mrisLDw+0NAAAEBlfNuFStWlXBwcHavXv3MR83j6tXr+7YuAAAgDu4KriEhYWpQ4cO+v7774s/ZjbnmseJiYmOjg0AADjPdUtF5lLowYMHq2PHjurcubOef/55paen69prr3V6aAAAwGGuCy4DBgzQ3r179cADD9gNuW3bttVXX331uw27AAAg8LjuHJeyug4cAAC4g8+e4wIAAPBnCC4AAMBnuG6PS0kUrXqZKScAAOAbil63T2T3il8Fl7S0NPsnR/8DAOCbr+Nmr0vAbM41Z77s2LFDUVFRKleu3Ck/T1FZo6kQYJNv2eB7Xrb4fpc9vudli++3b33PTRQxoaVmzZq2ozBgZlzMXzYhIaHUns984/mBL1t8z8sW3++yx/e8bPH99p3v+V/NtBRhcy4AAPAZBBcAAOAzCC7HYRqnH3zwQZqnyxDf87LF97vs8T0vW3y//fd77lebcwEAgH9jxgUAAPgMggsAAPAZBBcAAOAzCC4AAMBnEFyOMmvWLF144YX25D5z8u7HH3/s9JD82ujRo9WpUyd70nFcXJwuvvhirV692ulh+bWxY8eqdevWxQdEJSYmatq0aU4PK2A88cQT9nfLbbfd5vRQ/NZDDz1kv8dH35o2ber0sPza9u3bdfXVV6tKlSoqX768WrVqpfnz53vt6xFcjpKenq42bdro5ZdfdnooAWHmzJm68cYblZycrG+//VY5OTnq06eP/fcA7zAnS5sXzwULFthfLGeffbYuuugiLV++3Omh+b158+Zp3LhxNjjCu1q0aKGdO3cW33766Senh+S3Dh48qO7duys0NNS+CVqxYoWeffZZVa5c2Wtf06+O/C+pvn372hvKxldffXXM4wkTJtiZF/OiesYZZzg2Ln9mZhSP9thjj9lZGBMezS97eMfhw4c1cOBAjR8/Xo8++qjTw/F7ISEhql69utPDCAhPPvmk7Sd64403ij9Wv359r35NZlzgGikpKfbP2NhYp4cSEPLy8jRp0iQ7w2WWjOA9ZmbxggsuUO/evZ0eSkBYu3atXfJv0KCBDYxbtmxxekh+69NPP1XHjh3Vv39/+8azXbt2NqB7EzMucE2zt1n3N1OOLVu2dHo4fm3p0qU2qGRmZqpixYr66KOP1Lx5c6eH5bdMOFy4cKFdKoL3denSxc7eNmnSxC4TPfzww+rRo4eWLVtm99OhdG3YsMHO2t5xxx2677777M/5LbfcorCwMA0ePFjeQHCBa96Rml8srEV7n/mFvnjxYjvD9f7779tfLma/EeGl9G3dulW33nqr3cMVERHh9HACwtHL/WY/kQkydevW1ZQpU3Tdddc5OjZ/fdPZsWNHPf744/axmXExv8tfeeUVrwUXlorguJtuukmff/65pk+fbjePwrvMO6GGDRuqQ4cO9sousyH9hRdecHpYfsns19qzZ4/at29v912YmwmJY8aMsffNch28q1KlSmrcuLHWrVvn9FD8Uo0aNX73pqdZs2ZeXZ5jxgWOMTVZN998s12qmDFjhtc3dOGP3zFlZWU5PQy/1KtXL7s0d7Rrr73WXp57zz33KDg42LGxBdLG6PXr12vQoEFOD8Uvde/e/XfHWKxZs8bOcnkLweU3P+BHp/KNGzfaKXWzWbROnTqOjs1fl4cmTpyoTz75xK4979q1y348JibGngWA0jdq1Cg7lW5+ntPS0uz334TGr7/+2umh+SXzc/3bPVsVKlSw512wl8s77rzzTnv1nHnh3LFjh20rNgHxyiuvdHpofun2229Xt27d7FLR5Zdfrrlz5+rVV1+1N68x7dDwmD59umnK/t1t8ODBTg/NLx3ve21ub7zxhtND81tDhw4tqFu3bkFYWFhBtWrVCnr16lXwzTffOD2sgHLmmWcW3HrrrU4Pw28NGDCgoEaNGvZnvFatWvbxunXrnB6WX/vss88KWrZsWRAeHl7QtGnTgldffdWrX6+c+R/vxSIAAIDSw+ZcAADgMwguAADAZxBcAACAzyC4AAAAn0FwAQAAPoPgAgAAfAbBBQAA+AyCCwAA8BkEFwCutGnTJpUrV87e2rZtW+LnK3ouU7oHwHcRXAC42nfffafvv/++xM+zc+dOPf/886UyJgDOIbgAcDVTSGhuJVW9enVb4AnAtxFcAHjd3r17bXAwDbJFZs+erbCwsJOeTRkyZIguvvhi+1zx8fF26eeRRx5Rbm6u7rrrLtvmnpCQoDfeeMMLfxMATgtxegAA/F+1atX0+uuv28DRp08fNWnSRIMGDdJNN92kXr16nfTz/fDDDzaczJo1Sz///LOuu+46G4TOOOMMzZkzR5MnT9YNN9ygc845x34eAP/BjAuAMnH++edr+PDhGjhwoEaMGKEKFSpo9OjRp/RcZlZlzJgxNgANHTrU/pmRkaH77rtPjRo10qhRo+xszk8//VTqfw8AzmLGBUCZeeaZZ9SyZUtNnTpVCxYsUHh4+Ck9T4sWLRQU9Ov7LrNkZJ63SHBwsN0Xs2fPnlIZNwD3YMYFQJlZv369duzYofz8fHu586kKDQ095rG5zPl4HzNfB4B/YcYFQJnIzs7W1VdfrQEDBtilnWHDhmnp0qWKi4tzemgAfAgzLgDKxL/+9S+lpKTYvSn33HOPGjdubPenAMDJILgA8LoZM2bYw9/efvttRUdH2/0p5v6PP/6osWPHOj08AD6EpSIAXtezZ0/l5OQc87F69erZGZiTNWHChOMGo98qyR4aAO5FcAHgat26dbNdReaclpKoWLGiPaQuIiKi1MYGoOwRXAC4kjk4bu3atfb+qV42fbTFixcXXyoNwHeVKygoKHB6EAAAACeCzbkAAMBnEFwAAIDPILgAAACfQXABAAA+g+ACAAB8BsEFAAD4DIILAADwGQQXAAAgX/H/zyFjsrbY1egAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "T = 10\n",
    "\n",
    "## Uniform mesh\n",
    "model = TwoTankModel()\n",
    "submesh_types = {\n",
    "    \"left domain\": pybamm.Uniform1DSubMesh,\n",
    "    \"right domain\": pybamm.Uniform1DSubMesh,\n",
    "}\n",
    "var_pts = {x_left: 20, x_right: 20}\n",
    "param = pybamm.ParameterValues({\"D_l\": 2, \"D_r\": 6})\n",
    "param.process_model(model)\n",
    "# solve\n",
    "t = np.linspace(0, T, 100)\n",
    "spatial_methods = {\n",
    "    \"left domain\": pybamm.FiniteVolume(),\n",
    "    \"right domain\": pybamm.FiniteVolume(),\n",
    "}\n",
    "\n",
    "\n",
    "def solve_model_with_uniform_mesh(L_min, L_mid, L_max):\n",
    "    geometry = {\n",
    "        \"left domain\": {\n",
    "            x_left: {\"min\": pybamm.Scalar(L_min), \"max\": pybamm.Scalar(L_mid)}\n",
    "        },\n",
    "        \"right domain\": {\n",
    "            x_right: {\"min\": pybamm.Scalar(L_mid), \"max\": pybamm.Scalar(L_max)}\n",
    "        },\n",
    "    }\n",
    "    # mesh and discretise\n",
    "    param.process_geometry(geometry)\n",
    "    mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n",
    "\n",
    "    disc = pybamm.Discretisation(mesh, spatial_methods)\n",
    "    built_model = disc.process_model(model, inplace=False)\n",
    "    solver = pybamm.IDAKLUSolver()\n",
    "    t_eval = [0, T]\n",
    "    solution = solver.solve(built_model, t_eval)\n",
    "    return solution\n",
    "\n",
    "\n",
    "solution = solve_model_with_uniform_mesh(1, 3, 6)\n",
    "ax.plot(solution[\"x\"](t=T), solution[\"phi\"](t=T), label=\"Uniform mesh\")\n",
    "ax.set_xlabel(\"x [m]\")\n",
    "ax.set_ylabel(r\"$\\phi$\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Symbolic mesh\n",
    "\n",
    "We now show how to use the symbolic mesh to optimize the geometric parameters. The only difference in the two codes is that we must use the pybamm.SymbolicUniform1DSubMesh as the submesh type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x3463d7e10>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGwCAYAAACOzu5xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVupJREFUeJzt3QdYlWXjBvD7HPaWISgyHCCIC3EA7twjU8u0NHOXlrmqz2zY1raj3KW23LPMkZqmKaCgKDjBBYoyVLbM8/6v9yH9a1k5gOeM+3dd53vPOeLhhvg8N8/7vM+jURRFAREREZEB0MoOQERERHSvWFyIiIjIYLC4EBERkcFgcSEiIiKDweJCREREBoPFhYiIiAwGiwsREREZDHMYEZ1Oh5SUFDg4OECj0ciOQ0RERPdAXVIuJycHnp6e0Gq1plNc1NLi7e0tOwYRERE9gOTkZHh5eZlOcVFHWm5+4Y6OjrLjEBER0T3Izs4WAw8338dNprjcPD2klhYWFyIiIsNyL9M8ODmXiIiIDAaLCxERERkMFhciIiIyGEY1x4WIiB58OYmioiLZMchIWVhYwMzMrFxei8WFiMjEqYXl3LlzorwQVZQqVaqgWrVqD73OGosLEZGJL/x1+fJl8duwejnqfy3+RfQgP2P5+flIS0sTj6tXr46HweJCRGTCSkpKxJuKumKpra2t7DhkpGxsbMRRLS/u7u4PddqI1ZqIyISVlpaKo6WlpewoZORs/yzGxcXFD/U6LC5ERMT93chgfsZYXIiIiMhgsLgQERGRwWBxISIik1azZk3MnDnz1uMrV66gc+fOsLOzE5fwmvKpnQ0bNkDfsLiQyVN0OmRl56KguGySIhHpv/bt22PChAl/e37p0qX3XTYOHjyI55577tbjGTNmiEvEY2Njcfr06XLJS+WHl0OTSSgsLkba+eO4fjEBN9LOQLl2Dpa5F+FUcAnuJZexV9cYLysTEFrbFe3rVkVHr1L4+NSGhmtaEBm9qlWr3vH4zJkzaNq0Kfz9/R9qUT9eqVUx+K8yGa2SUh2WH0jCI5/tRtDUrXD/oQMa/T4CoSemISx1OULy9qJO6Vk4aG7AR5OKwhId9pxOx8ebYuG+JByX3w9E1FfDELtzBfJzs2R/OUSVt1hYUYmUm/q5y9vQoUPRp08ffPbZZ2LhM1dXV7z44ot3XJJ7+6ki9f7atWvx3XffiVMl6t9XJSUloXfv3rC3t4ejoyP69++P1NTUW6/xzjvvIDg4GF9//TVq1aoFa2tr8bz6GgsWLMCjjz4qLgeuV68eIiIikJiYKEaN1NNRLVu2FGXpn5w/f168zqpVq9CmTRuxJkrz5s3FaJA6WtSsWTORq3v37khPT7/j76p51M+p5gkMDMTcuXPvKFdjx44V3xf1z319fTF9+vQ7/n5GRgb69u0rsqtF7qeffoJsHHEhozz1c2TXKiREbcGU7P5/PqtFIrxho1WQZe2JQntvwLkmbDxqw8mzLvy962J7roLdp9KRFL8XZqml8FRS4ZmxDti7DoV7LBBn3RD5dbqjad8JMLfgb1JknG4UlyJo6jYpn/v4e11ha1n+b0u7du0Sb87qUS0MAwYMECVj1KhRf/tYtQg8++yzopzMmjVLlAR1K4SbpeX3338Xi/ap5Ud9nd27d9/6u+prq6Vn3bp1dyyw9v777+OLL74Qt8mTJ2PgwIGoXbs2pkyZAh8fHwwfPlwUiC1btvzr1/H222+LgnXz76iv4+DgIHKqxUItU1OnTsW8efPEx//444/i8VdffYUmTZrg8OHD4mtWy9KQIUMwe/ZsUUTUQqS+ZnJysrjd7t1338Unn3yCTz/9FF9++SUGDRqECxcuwMXFBbKwuJBRSTi8B0Vb3kBw0VEEA/jFJgRtOj6Gxxp7ws2+x7+uI+BvB/h7OABtayMvpzeOR25G4clt8Lm6D9U16WhYeAg4fgjrzp1CyIjZqOlmV6lfGxE9GGdnZ/HmrZYJddShZ8+e2Llz512Li3rayMrKShQWdV8d1fbt2xEXFyf2c1K3RVCpIzL169cXRUcd/bg5gqE+/9dTT8OGDROlQqUWl/DwcLz11lvo2rWreG78+PHiY/7LK6+8csffefrpp8XX0apVK/HciBEjxByf24vO559/jscff1w8VkeCjh8/LkaA1OKijiKpoyitW7cW/zaqIy5/pY44qZ9HNW3aNFF2Dhw4gG7dukEWFhcyCinnTiJl3etolrNTPC5ULHDIcwBm9R8EJ2e3+349O4cqCO48EOg8UIzgXDgdi8sH1sPj7GpMu94R+bP3YuqjQRjQ3JsLd5FRsbEwEyMfsj53RVALxu0jIOroi1pE7tWJEydEYblZWlRBQUFiErD6ZzeLi/rG/9fSomrUqNGt+x4eHuLYsGHDO54rKChAdna2GOn5J/fyOml/7geUl5cnTj+pZeb2gqaOFjk5Od0qJerVUwEBAaKIqKezunTp8o+fUx2pUfPd/ByysLiQQcu8fg0nV7yOkCur4akpgU7RIKZKZ3g/MQ3hPg8+se526gRd38AQcbt0bTL81sQh8uw1vLYuDqVRC9F9wBi4uNcol89FJJtaxCvidE15U99As7L+PvcsMzPz1hvzTRYWFn/7GitiJ2z1jf1ubv/8N3/Rudtz/5XpXl5H9+dr5ObmiuOiRYsQGhp6x+vcLHEhISFiFEk9RbVjxw4xKtSpUyesWbPmrp/zr59DFk7OJYO1Je4yOs/aB+8rv8JSU4I4qyY4+/gvaD5xNaqVU2n5qxou9lg2Mgyv9wjEY+aRGHT1S+jmtsSR31ZVyOcjortTRwkOHTr0t+fV5+rWrVuun0ud3PrX+R/qKRe1JKkjL/rIw8NDbJx59uxZ+Pn53XFTTxndXgDVuTpqwVm5cqWYo3Pt2jXoM/2v1UR3se7QRbyy+gh0ihYLXF/AE81roVHbvpVy+bJWq8FzbevgjHNPnF+/ATV1SXDbMwpRxzaj0fAvYWPnUOEZiEzdmDFjxLyVcePGYeTIkWJeyi+//ILly5fj559/LtfPpY5CqKdk1Imp6uRY9XTLCy+8gHbt2okrevTVu+++K74/6giUeiqosLAQ0dHRuH79OiZNmiQmC6unzdSJu1qtFqtXrxbzevR90T2OuJDBiVo7AzvWLoROAfo388LbL7+Mxu2fqPQ1V+o0DEO1VyIQ6V426S706nqkfx4qJggTUcVSr8rZs2cPTp48KYqFejpEvTpGffMt74mj6umRjRs3ikm+bdu2FZ9P/fzqCIU+GzlypLgcesmSJaJ4qUVLnbx7c8RFvSJJvWJILV/qPB31suvNmzeLEqPPNEpFXDgviTqxSW2W6nnPf5vgRIYrasV0hJ78CMWKGRYELsYLA3qLERDZ4vash8dvk+COayLbkRYfo1nPv1+xQKRv1Emh6jyH29ceIarsn7X7ef/W71pFdJvI76eK0qKKqdYfLw54TC9Ki6ph276weikSh+zaQAMFqyISsS8xQ3YsIiKjw+JCek+9HDliyWSEnZklHkd4DUfo83P1bjl+J1cPBE/aiM9852JVSVs8/30MjqVwxV0iovKkX//yE92ltER+PQHhF+aLxxE1xyB85Ay9Ky03ac3MMPHZ/gir7YLcwhJMWrwdKedPyY5FRGQ09PNff6I/90xZ/cN8hKd8Kx5H+k9C+NCyU0X6zMrcDAsGN8MjVXOxsGgKSr7ri8yMK7JjEREZBRYX0ks6nYI3N8Rj8nFvLCt5BFH1piBs0NswFE42Fvj4qRaw0urgo7uE1Pm9cSMvR3YsIiKDx+JCemnqT/H4MSpJXbYW5o/NRuiA12Bo3GvUQuGAVciCHQJKTuLknP4oKS6SHYuIyKCxuJDeif5lEYKi34KNphAz+gejfwsfGCrfek2R0n2J2DupSf5+xMwbIebtEBHRg2FxIb3bLDHgwFsYaL4LcwOOok8Tw98DqF5oVxxv+YXYRyn02k+IXGp4o0dERPqCxYX0hnoaJfvHoXDQ3MBJiyC0GfQ6jEWTrs/iYNAUcd/9/E9Yvf+k7EhE9AB2794tVtJV9yl6GDVr1hTbB9ykvuaGDRtgSGRlZnEhvXHw28kILDmBbNjCcdBSmFtYwpiEDpiMbbWnoF/R23ht01nsOil3a3giQ5aeni72K/Lx8RH7FKl77HTt2hX79u2DIbp8+TK6d+8uO4ZBYHEhvXBs/2aEJi8R9xOavw/PmgEwRl0GT0bHpkEo1Sl4efURXM0tlB2JyCA98cQTOHz4ML799lucPn0aP/30E9q3b4+rV6/CEKnFSy1g9N9YXEi6zIxUuP06FlqNggNVeqBpz5EwVurQ6rS+DVHPww49Cn7Bya9HyI5EZHDU0zR79+7Fxx9/jEceeQS+vr5o0aIFpkyZgscee0x8zPDhw/Hoo4/e8feKi4vh7u6Ob775RjxWi85LL72ECRMmiA0UPTw8sGjRIuTl5WHYsGFiE0I/Pz9s2bLlbxnUkZ1GjRqJPXfCwsIQHx9/x5+vXbsW9evXF2VEPS30+eef39dpl4sXL+Lpp5+Gi4sL7OzsxEaIUVFRd/276uaI6t9XN5ls06YNbGxsxKaJaqE7ePCg+Lv29vZiREcdqbqdugljvXr1xNcRGBiIuXPn3vqzoqIijB07Vuwgrf65+n2ePn36HX8/IyMDffv2ha2tLfz9/UWBrHCKEcnKylI3jBRHMgw6nU6ZtuBb5epUTyXpnUAlN/u6YgpOxkYqJVOdFOVtRyV2+zLZcciE3bhxQzl+/Lg43qEw959vRffzsfn39rH3obi4WLG3t1cmTJigFBQU3PVj9u3bp5iZmSkpKSm3nlu3bp1iZ2en5OTkiMft2rVTHBwclPfff185ffq0OKp/p3v37srChQvFc2PGjFFcXV2VvLw88Xd27dol3mfq1aun/Prrr8rRo0eVRx99VKlZs6ZSVFQkPiY6OlrRarXKe++9p5w6dUpZsmSJYmNjI443+fr6KjNmzLj1WH3N9evXi/tqvtq1aytt2rRR9u7dqyQkJCgrV65U9u/ff9ev9dy5c+LvBwYGKlu3bhX/PcPCwpSmTZsq7du3V/744w/l0KFDip+fnzJ69Ohbf++HH35Qqlevrqxdu1Y5e/asOLq4uChLly4Vf/7pp58q3t7eyp49e5Tz58+LLMuW/f+/V+rn9PLyEs+pGceNGyf+u1y9evX+ftbu8/2bxYWk+j7ivOI7eZPS6vXvldNHoxRTsn/eGFFcUt+uqWRdz5Adh0zUP76ZvO34z7cf+t35sR9U++ePXdzjzo/9uNbdP+4+rVmzRnF2dlasra2Vli1bKlOmTFGOHDlyx8cEBQUpH3/88a3HvXr1UoYOHXrrsVpcWrdufetxSUmJKDaDBw++9dzly5fF+0pERMQdxWXFihW3PkZ9o1aLiVouVAMHDlQ6d+58R5ZXX31V5LmX4rJgwQJRqP6pAPxTcfn6669vPbd8+XLx3M6dO289N336dCUgIODW4zp16txRRFRqeQsPDxf3X3rpJaVDhw7iF8y7UV//zTffvPU4NzdXPLdly5YKLS48VUTSnL6Sjfc3HRf3h3RtCf+GLWBKggd/jIua6nDHNZz4boLsOEQGN8clJSVFnJro1q2buNonJCQES5cuvfUxI0eOxJIlZXPnUlNTxSkf9RTS7dTTPTeZmZnB1dUVDRs2vPWcevpIlZZ252T68PDwW/fV0zkBAQE4ceKEeKweW7VqdcfHq48TEhJQWlr6n19bbGwsmjRpIl73fjS67Wu5mfuvX8vNr0M9HXbmzBmMGDFCnEa6efvggw/E86qhQ4eKLOrXNm7cOPz666//+jnVU1qOjo5/+16VN/MKfXWif1BwIw83vu6B9rpOyPfvjhGta8HU2Ng5IKvzF/D69Wmxvkv8vp/RoFUv2bGIyrye8s9/pjG78/Grif/ysX/5/XhCHMqLOu+ic+fO4vbWW2+JovL222+LN1zVs88+i9deew0RERHYv38/atWqJeaA3M7CwuLOuBrNHc+pj1W6Slw4Up2j8iAs7pL7r8/d/Dpyc3PFUZ3TExoaesfrqAVOpRbBc+fOicK3Y8cO9O/fH506dcKaNWvu+jn/+jkqCosLSXHkm5cQWhKHaZZJ0PUZC6227P9kpqZ+yx6IiumD0Ksb4LzjZeQHt4OtnaPsWESApZ38j71PQUFBd0xwVUdP+vTpI0Zd1PKiTrgtL5GRkeJSbNX169fFRFh1kqtKPf71smz1cd26dW+Vgn+jjmKok2avXbt236Mu90odffH09MTZs2cxaNCgf/w4dQRlwIAB4tavXz8xulWRue4FiwtVutgdyxGasVbcv9TuCzRydYMpC3p2BlJn7IWHLh1r1q/EU8+Mkh2JSK+plzw/+eST4rSP+iavXv0THR2NTz75BL17977jY9VRGPXqIvUUzZAhQ8otw3vvvSeKkVoA3njjDbi5uYmSpHr55ZfFVT3vv/++eMNXS9NXX311xxU7/0a9mmjatGni9dSreNSretRLv9WiEX7bKaqH9e6774pTQE5OTqKQFBYWiu+jWsQmTZqEL774Qnxu9bSVVqvF6tWrxWXbVapUgUwsLlSpsq6mwveP/4n7kR5PIeyRfjB1Dk4uSOwwCyO2nMfxY54ISLqOJj7OsmMR6S11LoZ6emPGjBliPoZ6mbO3tzdGjRqF11+/c8Vt9dSG+uarXpqsvvGXl48++gjjx48X81aCg4Px888/w9LS8tYpFvXS5KlTp4ryon5+tejcPIX1X9TXUeeTqAWoR48eKCkpEaNJc+bMQXlSS516GfOnn36KV199VcxRUefEqJeHq9RCqJZB9WtUR4rUMrZ582ZRYmTS/Dkz2ChkZ2eL5piVlSWGt0j/RM59HmFpK3Be64Pq/4uClbWt7Eh6Y+LKWKw/fAn+7vbYNK41rMz/e0iZ6GEVFBSIeQzq/A91zoixUedy1KhRQ5wuevzxx2XHMWkF//Kzdj/v37yqiCrNxcR4hKSuFvdz2r3L0vIXUx8Ngpu9JTTpJ7B7xQzZcYgMmjpBVL26RR3xUE9t3FyYjgyf3hYXdRhOnZ18c8iKDF/k1h9gqSnFUetmaNiOv/n8lbOdJT57xBabLF9Hh4RpOBt/91Uyiei/JSUlifkny5Ytw+LFi2FuzpkRxkIv/0uqSxQvWLDgjuvDybDFXLiGVy62wUqtOz4e0F52HL3VLjwch/eFIyT/D5SufwElARFGt9kkUWVQl9k3opkQpM8jLur5SPXSLPXacnXvCDJ86j8eH/xStjBTnaadULteiOxIekuj1cLnmbnIhh38SxMRveID2ZGIiPSK3hWXF198ET179hQzwf+LeumWOqHn9hvpn/2/b8PlpDOwsTDDpM51ZcfRe26evjjZ+DVxv3HiPGRcSZIdiUwARyfIUH7G9Kq4rFixAocOHfrb7pP/RP04dRbyzZt6ORzpl8KCfNT8fRx2W03C+42vwt3R+K5aqAjNe4/FKfNA2GiKkLj2PdlxyIjdXBBN3QmYqCLl5+ffdbVdg53jkpycLK6J3759+z1fkqduYa4uknOTOuLC8qJfDq/9DGFKKtI1zujR7c4t5unfTxkVtX0N+G0omqRtQOrFKfDwqiM7FhkhddKqupZHenq6eEORvUYHGedIS35+vrjKS73C615WDzaI4hITEyO+KHXhnpvUlQ737NkjVhxUTwv99Yu1srISN9LfxebqJcwX9881nIAW9k6yIxmUBq17I2Jfa+zI9UVJ1HW86yU7ERkj9epNdYE0dX2NCxcuyI5DRqxKlSpi5d2HpTfFpWPHjoiLu3PzLXVficDAQEyePPmhGxpVvhOrpiIMeTinrYmmvcfKjmOQoy6aAd/hm4WRsDiUhpGP5MPbhWvfUPlTV2r19/fn6SKqMOpoXnm9j+tNcVGXFm7QoMEdz6nLD6t7Qfz1edJ/l84eQ8iV1YAGyG33Nsy4hsIDCavtilZ+rtiXeBVf7jyNT54Mlh2JjJR6isgYV84l48OTmVQhUtdN4WJz5WRS5wB01R7AkLhnxerDRESmTK9/Dd69e7fsCPQAYs5nIDLTCYFmVnDoNU12HIPX1NcZLzruQ/2CC4j+6R14TVojOxIRkTQccaHyX2xu8yl8WvIUPqm/HrXqh8qOZBRsu04Vx5CsHbhw8pDsOERE0rC4ULnaHHcFh5MyxWJzL3RrKjuO0fBr0haxti2h1SjI2PSu7DhERNKwuFC5LjZntXEUQjSn8Xy72lxsrpw59nhHHJvm7uYGjERkslhcqNzE/jwXnUr3Yr7VbDzXiouOlLfaDUIRY1+2QWXWZo66EJFpYnGhcqErLYXn8W/E/bN1R8DWhuuNVAS3R99GqaJBk/x9OB13QHYcIqJKp9dXFZHhOLJzOZooKWJX44a9uNhcRfENDMHmqsOwKsUVmoMKljSUnYiIqHJxxIXKhdXBeeJ4zPMJ2DlUkR3HqAU99QH2IgS7Tmcg5sJ12XGIiCoViws9tNOHdiOoOB5Fihn8Hn1ZdhyjV9PNDk+E1BD352yLlR2HiKhSsbjQQ8v5bYY4HqnSGVU9a8qOYxJe6uCPMRabMOPSQByL2CI7DhFRpWFxoYeSfC0fq6/744yuOlw7T5Idx2Somy22d8+HkyYfmt8+gKLTyY5ERFQpWFzooSzedw4rSh/BOz5LxeW6VHlq9p2KQsVCnKY7EbVNdhwiokrB4kIPLCu/GCsPJov7o9rWkR3H5Hh41cER127ifsG+ssnRRETGjsWFHljMT3PwWOl2NPSwQht/N9lxTJJbh5fEsVHOXqRdOic7DhFRhWNxoQdSVFiABidn4SOLr/FWrdPQaDSyI5kk9fTcccuGMNfocHbLl7LjEBFVOBYXeiBHtnwDd1xDOpwR3H247Dgm7UbwCHH0u7gWhUWFsuMQEVUoFhe6b+oVLK5HF4r7ibUGwdKKmynK1KjTQHyr7YsBhW9g6/EM2XGIiCoUiwvdt/g/NqK27jzyFSsE9ZogO47Js7C0QlarN3BGqYFv95+XHYeIqEKxuNB9U/Z/JY5H3R+Dk0tV2XEIwFMtvGFhpsGhpEzEJWfKjkNEVGFYXOi+nD1+AI0KosUOxT49XpEdh/7k7mCNYXUL8aXFbOSt4yaXRGS8uDs03Zd1MRfRtLQxnJyqIKRWoOw4dJvHgxwQeC4ShdcskJlxBVXcqsmORERU7jjiQvcsNbsAC05aY1jxZGj6fS07Dv1FQNMOSDSrAytNMU5uniM7DhFRhWBxoXumTvwsLlXQvKYzmtR0lx2H/kKj1eJag6Hivu/Z5SgtKZEdiYio3LG40D3Jy8lElciP4YFrGNmmtuw49A8adR2OTNijOtIRt2uV7DhEROWOxYXuSfwv8/Ac1mGFzcfoFMjRFn1lbWuPE9X7iPtm0YtkxyEiKncsLnRPC865n14h7qcFDIKZGX9s9Jlv13HQKRo0LDyEC6diZcchIipXvKqI/lPi0X3w151HoWKBwC4jZceh/+BZMwA7nXoj4qodzI8V4LUA2YmIiMoPf3Wm/3Ttj8XiGO/YhgvOGQjLXp/h69Ke+OFIFnILOUmXiIwHiwv9q4IbeaiXsU3ct2w2WHYcuket6rihdlU7UVrWH7ooOw4RUblhcaF/dey35XBEHq7ADUGtHpMdh+6RVqvB0Bae6K39A9V/GyfmKRERGQMWF/pXJxISkatY45zXYzAz55QoQ9K3kRumWSxGp+LfcWz/JtlxiIjKBYsL/aPLWTfwVmobtCicC6/u3JfI0Dg4uSC+ag9xvzhigew4RETlgsWF/tG6Q5egKECDWp7wrlFDdhx6ANU6lm242Ch3H64kJciOQ0T00Fhc6K7UORGxUbvVe3iyqZfsOPSAfOs1RbxVMMw0CpLXvik7DhHRQ2Nxobs6cWA7FhW8jPVW76JHA+4ybMjMO70lFqRrnrUVR37jNgBEZNhYXOiuciOXimOxsx/srC1kx6GHENi8Ew5Ue0rcd98zBVm5ebIjERE9MBYXuuuGig2u7xT3HcLLdhsmw9b42U+xz6w5Xih8CR9uTZQdh4jogbG40N8c2/E9bDWFSNZ4IrB5Z9lxqBzY2DnAcvAqxMIfq6IvYvepNNmRiIgeCIsL/Y398bINFS/59oVGyx8RY9G8pguGtawl7s9fswXZmVdlRyIium98V6I7XEyMR1BxPEoVDWp1GiE7DpWzV7sGYIzjfnxb9DJOfjtOdhwiovvG4kJ3OPf79+J4zKYpPLzqyI5D5czG0gy9OrSBBUrQ4vomxP2+TnYkIqL7wuJCt5TqFLyW2hFDiiYjJ3Si7DhUQYLCuuGAez9xv+quV5GTdU12JCKie8biQrfsS8zApexixFo1Q9PW3WXHoQrUaMjnuKTxQDVk4Pi342XHISK6ZywudMvq6GRx7B3sCWsLM9lxqALZ2jshs/NMcT/02k+I27NRdiQionvC4kJC1rV0/O/UALxmvhz9gz1kx6FKUL9lD0S5PSHuV/3tZeRmX5cdiYjoP7G4kHBy+2J4a9LQxTIO9b1dZcehStJgyBc4r6mB74o74OMd52XHISL6TywuJLgkrBbHNL9+XLvFhNg5VMHlgTsxt7Q3vj+QIuY5ERHpM75DEc4dPwj/kgQUK2ao22m47DhUycL9q2NwmK+4/8bqGOTmZsuORET0j1hcCKm/fyOO8fbhcHGvITsOSfBa90B0cErBghuTcGwprzIiIv3F4mLiSktK4Je6RdzXNHlGdhySxM7KHJPaeCBAexGhGesQv+9n2ZGIiO6KxcXEnYreATdkIht2CGrTV3YckqhB68cQ5dpH3HfbPhEp507KjkRE9DcsLiZuZ5KCJSVdEeP6GCytrGXHIcnqD5mJi5rqqIZ0mH/bHedPRMuORER0BxYXE6bTKfgx0RLvlgyBrtO7suOQHrB3dIbVqK04r/WBO66hysreOB2zS3YsIqJbWFxM2OHkTFzJLoC9lTla+7vJjkN6oqpnTTi/uAOnzANQBbnI+Okt7EtIlx2LiEhgcTFhZ/csR5j2ODoFusLKnEv80/9zcvWA1/hfsd2uF8YUvoRhS6OxNf6K7FhERCwupkrR6dD6zBdYYfkBnnXhJEy6++J0bSd+i5YN/FBUqsMLP8Zg287tsmMRkYljcTFRCbF/oDrSka9YIahN2ZUkRH+ljsR9+XQTDGjmjWe029B1bz9E/sj5UEQkD4uLibp6cJU4nnAIh7Wtvew4pMfMzbT46ImG6O5TKh6HJXyBiEUTxKgdEVFlY3ExQeobjveVP4f8gzjaQv9No9EgbNSXiKg1VjwOv7QEB+YMg66kRHY0IjIxLC4m6Gx8JLyUK7ihWCKQi87RPVI33wwf8iGi6r8FnaJB6NUNODzrSRQVFsiORkQmhMXFBKVF/XmayD5UTMAkuh+hT76Cwy0+E5tyNs35DSdmPIobBUWyYxGRiWBxMTGKosDySoy4r6v3mOw4ZKCa9hyJE+0XilG7zTl+GLwkGlk3imXHIiITwOJiYk6n5uKJvMl4ouRDBLR9UnYcMmCNHumHhCd3YplFX0RfuI6nFkYiPadQdiwiMnIsLiZmc9xldbYCnP1D4eDoLDsOGbhGDRph5XPhcLO3QvLlK4ib2Rcp50/JjkVERozFxcRsi7sojt0bVJcdhYxEkKcj1owOx6d2P6BD6T6YL+2GCyfKTkcSEZU3FhcTcuHUYSzPGoxpFovRKdBddhwyIjXd7NB0+Eyc13qLzRmdVj6GU4d+lx2LiIwQi4sJSdm/Es6aXNS3y4KTnaXsOGRk3GvUQpUXduC0eV2xOaPXxv6I/+Mn2bGIyMjoVXGZN28eGjVqBEdHR3ELDw/Hli1bZMcyGu4Xt4ljUd1esqOQkariVg2e435FvFUw7DQFqLt9GA7/+oPsWERkRPSquHh5eeGjjz5CTEwMoqOj0aFDB/Tu3RvHjh2THc3gXUyMR53SsyhRtPBvO0B2HDJi9o7O8J+4GYftWsNSUwKPfVOxNipBdiwiMhJ6VVx69eqFHj16wN/fH3Xr1sWHH34Ie3t7REZGyo5m8JL3rxTHE9bB4rdioopkZW2HhhPWY69LPwwt+h9eXn8aX+89KzsWERkBvSoutystLcWKFSuQl5cnThndTWFhIbKzs++40d25JZWdcsv3e1R2FDIR5haWaP3S12jbqq14/MEvJ7BkwzZuzkhExlVc4uLixCiLlZUVRo8ejfXr1yMoKOiuHzt9+nQ4OTndunl7e1d6XkOgrqvhX5KAUkUDP54mokrenPGNnvXwatcAhGpO4OnDg3Bg7gjoSst2miYiMvjiEhAQgNjYWERFRWHMmDEYMmQIjh8/ftePnTJlCrKysm7dkpOTKz2vIfgtIRNflvTBHrvOcPXwkh2HTLC8vPiIHyaEaGGJEoRmrMOhWf1RXMRVdono/mkUdfMaPdapUyfUqVMHCxYs+M+PVU8VqSMvaolRr0qiMn3n7sPhpEy817s+ng2vKTsOmbCYTYvQ6OBkWGhKccQmFHXHroWNnYPsWEQk2f28f+vdiMtf6XQ6MZeFHkxK5g1RWjQaoGt9TsoluZo+OgrH281HgWKBxjeicG5mN2RnXpUdi4gMiF4VF/XUz549e3D+/Hkx10V9vHv3bgwaNEh2NIN19Pd16KSNQbiPHTwcrWXHIULjDv1xrsePyIYtgorjkf5lJ2Skp8qORUQGwhx6JC0tDc8++ywuX74shozUxei2bduGzp07y45msHzj5+Bry2OIdFI7anvZcYiEeqFdccZmLYrXPYXjRe6YuTQe3450gJezrexoRKTn9H6Oy/3gHJc7ZVxJhsu8htBqFFwZEYNq3n6yIxHdIfnsSTyz8gIuZJWgmqM1fhjZAn7unPNCZGqyjWmOCz24c5EbRWlJNKvD0kJ6ybt2IFa80BZ+7vZIzc7HwbmjkHB4j+xYRKTHWFyMmDZxuzhmVG8nOwrRP6ruZINVz4fjLZddeBpb4LnhScTv2yQ7FhHpKRYXI1VSXAT/3APifpXGPWXHIfpXLnaW6D/6LRyzbCw2Z/T/dSgOb/9Rdiwi0kMsLkbqdMxvcEQ+MmEP/yaclEuGsTljnYlbEGvbElaaYjT8YywObpgjOxYR6RkWFyOVfux3cUx0DIOZuV5dPEb0j6xt7NBg4kYcrNIN5hodmse+jshlH8iORUR6hMXFSH2U2wPtCz9HZrPxsqMQ3ffmjE1fWoZI97J9tZqcmoGvf94FI7oAkogeAn8VN0JXsgpw4nI2NJrqCGkaJjsO0X3TmpkhdPR8RH7rjCUJ1ti2Lx/JumN4u1d9aLUa2fGISCKOuBih3afSxLGxVxUx6ZHIEGm0WoQN+xitew0TW1Z8G3EB7/24jZszEpk4Fhcj5LP3FcyzmIH+Na7JjkL00AaH18TMAcHw0WbgucQxODbjMRTk58qORUSSsLgYmaLCAjTK3oPuZgfR3JerB5Nx6B1cA7O6OMAFOQi+EYmz3JyRyGSxuBiZ09HbYa+5gWtwRJ1GrWXHISo3Tdo/jnPdvkeOYoOgojikfdkZV1Mvyo5FRJWMxcXI5MZtEcczTuFigiORMakX3h2pj68Vxdyv9Azy53fGlaQE2bGIqBKxuBiZamll+7xoA7rIjkJUIfwat0LeoE24gqrwVlKgWdwVFxLiZMciokrC4mJELl84hZq6ZJQqGviF9ZYdh6jCePs3BkZsxQWtF67q7DF4+RnEXcySHYuIKgGLixFJitoojqctg+DkUlV2HKIKpe547jhmB6a7foCkfAs8vSgSEWc4YZfI2LG4GJFDGVpE6+riuldH2VGIKoVz1eqY+3x3hNV2QW5hCbYvfQ+HdyyXHYuIKhCLi5EoKC7F7MsN0K/oHVTp9IrsOESVxsHaAkuHtcBLNS9iqtlSNNz7Ag5unCc7FhFVEBYXI3Hg3DXcKC6Fh6MV6nly/RYyLdYWZhg/YhgOOnUp25zx8GuIXD5NdiwiqgAsLkbixOF9cEQuHglwh0ZdH53IFDdnHLcCkVWfFI/DTn2MiMWvQtHpZEcjonLE4mIkepycgkNWo/G48xnZUYjkbs44ZiEifJ4Xj8OTFiJq3nPQlZbKjkZE5YTFxQhcTIwX61ko0KBe03ay4xBJ35wxfPgniAyYLB6Hpa/G10sXobiUIy9ExoDFxQhcPPiTOJ62agAHJxfZcYj0QtjTryM65CPMLn0c0xK8MeaHQ2ISOxEZNhYXI2B74TdxzPF5RHYUIr3S7LExCHr6I1iZa7HjRCrGfLMLOVncNZ3IkLG4GLgbeTkIuBEr7ldv2kt2HCK90ynIA98ObwE3q1KMTnkDV2Z3xrW0S7JjEdEDYnExcAlRm2GlKRb7tvgEhMiOQ6SXwmq7YvmTnqirTYF/aSJy5nfBleRE2bGI6AGwuBi4G8e3iuMF11ZiUiIR3Z1/g2bIfvpnXIEbfHUXgW+6Ivl02WglERkOvtMZMEVR8FFud0wpHgFN8CDZcYj0nm9AMDB8K5I1nqiGDNgt64XEI3/IjkVE94HFxYCdSc/D4Uw7rEVnNAjtIDsOkUGo5uMP29HbkWhWBy7Ihse6fjgWUTZySUT6j8XFgO0+lSaOobVdYGtpLjsOkcFw9fCCx7gdOGbZEMUww6u/JOO3k6myYxHRPWBxMWBVD3yCZ822oVstS9lRiAyOuuZRnQlb8UWNWThe4onnvovBhsO82ohI37G4GKi8nEx0z16F9yy+RVtv/mckehDWtvZ4e8Tj6NukBkp0Clau/hFRKz+WHYuI/gXPLxiohMjNCNaU4qKmGrzqNJQdh8hgWZhp8fmTjeGrzcDI+M9hf6IAEUuuIWzIdF6pR6SH+P9KA1Vwumy13EsuYfzHleghabUajH+iA+J8BovH4RfmI2r+aG7OSKSH+I5noDyuHhBHC7/2sqMQGc/mjCM+Q2TdV8XjsLSViJk9ECXFRbKjEdFtWFwM0NXUZNTSXRD3azXrJjsOkVEJG/gmDgZPQ4miRfOsrYib0RsFN/JkxyKiP7G4GKDzMdvE8YxZLThXrS47DpHRad7nRcS3noNCxQJN8vdj3dw3kVtYIjsWEbG4GKa05EQUK2ZIdwuVHYXIaAV3HojELt/iF6UVpqa3x8BFkbiWx9NGRLKxuBigT3K6oVHhIhSFT5Adhcio1W/VE96jlsHB1gZHL2bhqfl/IPVykuxYRCaNxcXApGTewLmMPBRpbdAk0E92HCKj18irClaPbonqjlYYfH0OlAUdkJxwRHYsIpPF4mJgIhLTxbFhDSc4WFvIjkNkEvzc7bFmWBDaWRxHNaTD9sdeOHN0v+xYRCaJxcXAeP8xGT9bvo5BLqdkRyEyKTWqe8Luz80ZXZEF97WP43gkN2ckqmwsLgZE0engk3kADbXnUa+Gi+w4RCa5OaP7S9tx3KIhHDQ3UHvLMzjy2yrZsYhMCouLAbl0Vh2mzkCRYoY6TTvKjkNkkhyruKLWhK04YhMKa00xgn4fjehNC2XHIjIZLC4G5NLhsvVbEq2CYGPnIDsOkcmysbNH0MSfEe3YCRooWBhxGd9Hli0KSUQVi8XFgJgn/SGOWdXCZUchMnkWllYIGb8KiwPm4dfSZnhrQzy++i0BiqLIjkZk1FhcDGh+S82cGHHfKYiniYj0gdbMDKOeHoBxHcqWJlix/Q/sWfSK+P8rEVUM8wp6XSpnF04dQk1k4YZiCb8m3FiRSF9oNBpM6hIAFysF7XZORK2UVByYnYKQF7+FuYWl7HhERuehRlyOHj2KOXPmYNGiRTh27Fj5paK/iU26jg2lLRFj1waWVtay4xDRXwxtF4iMkHFic8YWmZsRN7MvN2ck0qfiMmvWLAQHB+ONN97Aa6+9hoYNG6Jx48aIjY0t34QkbE51xoTisYhr8ansKET0D5r3GYu4Vl+iSDFHk7w/kDizB3Kzr8uORWS6xWXx4sU4dOgQCgsL8eGHH+Kjjz7C9evXcfXqVZw9exbdu3dHmzZtsH8/V5QsT6U6BZFnr4r7Leu4yo5DRP+iSZdncLrzEuQp1mhQGIuU2V2QmXFFdiwi0ywun332GUJDQ2Fvby/KysGDB8XIy++//w5nZ2dRZNTbK6+8UnGJTdDpxER4FJ6Hg5UZ6ns6yo5DRP+hQevHcKn3KlyHA+qWnEbM/OdwJatAdiwi0ysux48fR05OjhhRsbCwgFarxYoVK9CjRw+4uLigdu3aWL9+PWJiYvDLL7/g/PnzFZfchGRFfY/tVv/DIvsFMDfjhWBEhqBuSDtkP/UTDmoa4X+5T+GJefvFBqlE9HDu+13Q2toazZs3R6tWrcSclsjISFFm4uLi8MEHH8DPzw/FxcV49tlnRZFxdOQIwcOyvfTnqbcaIbKjENF98A0MQfVx2+Do5olLmTfw5PwInDrLheqIHsYD//r++eef45NPPsHIkSPFvBd/f3/06tVLnEby9PQUp5KSkpKwahX38XgYxUUF8LtxVNx3b9RZdhwiuk9ezrZY9Xw4gqo7ouONraj+XThORv0qOxaR6a3jol5RpJ4SGj16NMLCwm6tFmlubi4m8aq8vLzEjR7cmdg9CNQUinPlNes1lx2HiB5AVQcrLB/VAskzX4djUR4sNg/Ckfx5aPxIP9nRiExrAbo6depg+/btSE1NFaeMioqKEB4ezrJSjjKP7RTHc/YhCDEzkx2HiB6Qk60VLMdvxpEv+6JxwUEE7X4O0fnX0aznKNnRiExv5VwPDw/07t27PF6K/sLhcoQ4Fvu0lh2FiB6SujlqvYmbEP3VQDTL2YmQA68iKi8Tof1flR2NyGDwEhU9VpCfC7/C4+K+Z3AX2XGIqByoK1+HTFiNKNe+0GoUhB7/AJFLX+fmjET3iMVFjx2+mINRxZOwyGwAvPwayY5DROW4OWOLFxcjwmu4eByRcAXTNp9geSG6B9xkUY/tO5eNPbrGcK3bAxotOyaRMVH/Px0+cgZ+2RiGWRF2wN5zyLpRjGl9G3K9JqJ/wf936LH9ZzLEMZzL/BMZrZ69B+DTfo2h1QA/Rydi61fjUViQLzsWkd7iiIueUjdm6355Lhy1QWhZu73sOERUgZ5s5g1HGwuYrRyITtdjEDcjDrXHboCdQxXZ0Yj0Dkdc9NSZ6B0YZbYJH1p9Cy8XO9lxiKiCda1fDdU6jxebMzYsPIxLszpzc0aiu2Bx0VP5p38Tx0tVuOgckalo0Kb3HZszZs7tjPQU7vlGdDsWFz1VNT1KHLV12smOQkSVvTnjgI1Igwtq6pJQvKgzLibGy45FpDdYXPRQ1tVU1C45K+7XbNZNdhwiqmS+9ZqiZOhWXNRUh6eShhs/PI0TKZmyYxHpBb0qLtOnTxc7Tzs4OMDd3R19+vTBqVOnYGrORP8qFqY6r/WGWzUf2XGISALPmgGwfv5XxJnXx4TC5zFgYRRiLlyTHYtIOr0qLr///jtefPFFse+RugdScXExunTpgry8PJiS4sTd4pjqwvktRKZM/cXFZ9Ju2PiEILugBIO+jsIf8WWjsUSmSq8uh966desdj5cuXSpGXtRdqNu2bQtTYX29bJTJonYb2VGISDInW0t8PyIUY36MQfbpfWi4egRikt5H0x4jZEcjkkKvRlz+KisrSxxdXFzu+ueFhYXIzs6+42bocgqK0TfvNXQo/AxezXvJjkNEesDG0gwLBzfDJPcYOGny0CTqZUSt/lx2LCIp9La46HQ6TJgwAa1atUKDBg3+cU6Mk5PTrZu3tzcM3eGkTOgUDYqd68C9alXZcYhIT1iaaxH+0lJEufYp25zx2HuI+O4t2bGIKp3eFhd1rkt8fDxWrFjxjx8zZcoUMSpz85acnAxDF33hujg28737KBMRmS4zc3O0eHEJImoMFY/Dz85G5IIXoOh0sqMRmXZxGTt2LDZt2oRdu3bBy8vrHz/OysoKjo6Od9wMXYtDk/GVxWx0ck2XHYWI9HVzxlGzEOk3UTwOu/wjor98BqUlJbKjEZlecVG3dFdLy/r16/Hbb7+hVq1aMCXFxUUIyd+HR80iEeTJPUqI6J+FPfMODjZ6D6WKRmwNMGHFIRSWlMqORWRaVxWpp4eWLVuGjRs3irVcrlwp26dDnb9iY2MDY3f+WBT8NYXIhh18A0JkxyEiPdf88fE4UMUHE3dqkRufjsxvo7FgcFPYWurVP+1ExjviMm/ePDFXpX379qhevfqt28qVK2EKrh7/XRzP2dSH1sxMdhwiMgAtOvTF/KGtYGtphr0J6Vg961Wx+jaRsTLXt1NFpswi5YA43vDgwnNEdO9a+7vhx5GhiFj8PwzJW4Vzc3aieOTPcPP0lR2NyLhHXEyZelWAT+5Rcd8xgAvPEdH9aeLjjO79n0M6nFFLdwGFizrj0tkTsmMRlTsWFz2Rcv40quI6ihQz1G7M4kJE969WUHMUD9mKSxoP1FBSYfldd5w7VrbTPJGxYHHRE6fOn8cRXW0kWgbC2tZedhwiMlCetQJhNWo7zmlril+GXFb3xcnonbJjEZUbFhc9sT2zBnoXfYANwYtkRyEiA6fObXEZuwMnzevBCXnw/PkZ7D92RnYsonLB4qInbm5X37Smq+woRGQEnFyqwmfCNhyxbo43i4dhyLJT2Bx3WXYsoofG4qIHMrNzkJR6Vdxv5ussOw4RGQlbeyfUe3kbShv0Q3GpgrHLDmHN/pOyYxE9FBYXPXAh+hfEWY3EYvu5cLW3kh2HiIyIpYUZZj/VBE+38IGbch3Nt/ZCJDdnJAPG4qIH8hP3wUJTCgd7w99riYj0j5lWg2l9G+DDwLPw1aYhTN2cceFL3JyRDBKLix5wSj9UdscnTHYUIjJSGo0GnYdORWSd8eJxWMp3OPDVEG7OSAaHxUWywoJ81C46Je57NGgnOw4RGbmwwe/hQMN3xOaModd+wpGZj6OosEB2LKJ7xuIi2bm4/bDWFOM6HOHt10h2HCIyAS2emIgjYTPEgpchub/j5IyeyM/Nkh2L6J6wuEh2/eRecTxv2wAaLf9zEFHlCOk+DKc6fo18xQo2+Zcxeul+ZOUXy45F9J/4TimZ9eWyjRULPVvIjkJEJqZh28eR1GsFXtC+hT0XdRiwMAJp2TxtRPqNxUXybtgbbwRjU2koHOt1kh2HiExQYLMO+HJ0T1R1sMLJKzlYMOcjpJzjWi+kv1hcJDqbkYelN1rjZWUi6jRuKTsOEZmowGqOWDu6JQY6HsUbBTNh8W03nDt+UHYsorticZEo+nzZMv+NvarAytxMdhwiMmE+rraYNHQgLpj5lG3OuKo3Tkb/JjsW0d+wuEh09dhu1NJcRjPfKrKjEBGVbc744nacMg8UmzP6/PwU4vZslB2L6A4sLhL1ujAdu6xeRherONlRiIgEJ1cPeE/4FXFWIbDVFCJg53Ac2rpUdiyiW1hcJLmaehHeSoq4Xyv4EdlxiIju2Jyx7sRfcMi+HSw1JWgcMQG/7tguOxaRwOIiyYXYXeJ4Xusjtp8nItInVta2aDxhHQ4698Ti0u54bkchFu45IzsWEYuLLEVn94ljqnMT2VGIiO7KzNwczV76ARkt31R3O8K0zSfx+eYj3JyRpGJxkaTK1cPiqOXGikSkx9QVvaf0qI/J3QJhhSKERozBgTlDuTkjScPiIsGNvFzULk4Q92s06iA7DhHRfxrTvg7mtylES+1xhF7diNhZ/bg5I0nB4iLB2aN7YakpRRpcUN23ruw4RET35JGeT+Fw6Odic8amObtwcsaj3JyRKh2LiwT7cqtjRNHL2OwxmhsrEpFBadpjBE4+skhsztio4CCSZnVD1rV02bHIhPBdU4L9F4uwU9cUuob9ZUchIrpvjdo/gaSePyIbdggsPo6rX3VCxpUk2bHIRLC4VDKdTkHMhevifvOaLrLjEBE9kMAWnZHRbz0yUAVupWl49dudSL6WLzsWmQAWl0p2LiEOo0qW4RHLE6hX3UF2HCKiB1a7QSgKn92M/1m/hV3X3dFv/n6cTs2RHYuMHItLJcs4sg3jzDfgZZufYW7Gbz8RGbYatevj3ReHo66HPVKzCzFt/hKciilbYJOoIvCds5JpL0aJY457M9lRiIjKhYejNVY9H45e1bMwWzcd3j8NQNxebs5IFYPFpZLVyDkqjvZ+rWRHISIqN1VsLfHR8B64YB1QtjnjjuE4tO072bHICLG4VKJraZfgqaSK+z6N2sqOQ0RUruwcqqDuxM04ZNe2bHPG/eNwYN1s2bHIyLC4VKLk+LL9iS5oveDk7CY7DhFRhWzO2GjCWhys0hNmGgUtjr6FyB/elR2LjAiLSyXKP3dAHNMc6suOQkRUYcwtLNFs3A+IrDZIPA5L/AI//zgLiqLIjkZGgMWlEllmHBdHnWeI7ChERBVKXRU89LmvEFnrRewtbYCX43zx5oZ4lOpYXujhsLhUEvU3jVEFL6Fz4Sewb9JPdhwiokopL2FDpiG5x3co1ljgx6gkTFxxCEWFhbKjkQFjcakkSdfyce2GDhe0PvCvXVt2HCKiSjMwvA5mP9UE5lqg0fFPcWJGT9zI40J19GBYXCpJbHKmOAZ5OsJS/X8vEZEJ6dXYEz887o5BZjvRuOAgzs/sgqzrGbJjkQHiO2glsTv4JWZYzEHvKudkRyEikiKsWXNc6KFuzmiLesXHkSE2Z0yWHYsMDItLJfFK3Y2+ZvvQwCFXdhQiImkCQ7sg/YmyzRnrlJ5DwYLOSDl/SnYsMiAsLpWguKgQtYoTxf1qQVwxl4hMW52GYSh45hekaNzhpVyG+dJuuHAiRnYsMhAsLpXg/PGDsNIUIwt28KrNNVyIiLz8GsB85K84r/WGm3IdX67adGsuING/YXGpBFdPR4jjBet64vJAIiIC3GvUQpUXduATpylYc6MZBi2KxP5ETtilf8d30UpgllI2BJrn1lh2FCIivVLFrRrGvvgyWvm5Iq+oFJOXbEXMjhWyY5EeY3GpBFWz48XRplYL2VGIiPSOvZU5Fg9tjr6BtlhsNg3Be0fjwPqvZMciPcXiUsFy8vKRV6JFqaKBV4PWsuMQEeklK3MzfDqoFa67NC7bnPHIG4j48X3ZsUgPsbhUsLjL+ehZNB2drZfDzcNLdhwiIr3enLH5uB8R6fG0eBye8Bkivp4ERaeTHY30CItLBYu9WDZLvp6Ph+woRESGsTnj83MRWfNF8Tj84jeImjsCutJS2dFIT7C4VLCjSdfEsbG3k+woRESGsznj0GmICnoDOkWDsIx12D53PIpLOfJCLC4V7o2zg7HOcipCq3BDMSKi+xHa/3843PxTnFE88ealMDz3XTRuFHHkxdSxuFSgtEvn4I0raKw5A//atWTHISIyOE0fHYXk/tuRY+GKXafSMWTxAWTn35AdiyRicalAF+P/EMcLZr6wteepIiKiB9G+vhe+HxEKB2tzVEv6Gamft0ZG6kXZsUgSFpcKVHjhgDhmODWQHYWIyKA1r+mClcOD8YblcviXJiJ/fmdcvsDNGU0Ri0sFcrh6tOxOjaayoxARGbwgHw8UPfMTLqMqfJQUaJd0x4VTh2XHokrG4lJBSktLUbOg7LcBt0DuCE1EVB68/RpCO3IbLmi94YGrcFzeCwmxe2XHokrE4lJBLp6Ohb3mBvIVK/gENJEdh4jIaHh41YHjmO04bV4XzshB9fVPIn7fJtmxqJKwuFSQk5czsbW0OY7YtBCrQRIRUflxrlodnuN+RbxVsPglcc+2Nfj12BXZsagSsLhUkL3Z7hhdPBG7Gn0qOwoRkVGyd3SG34TN+MFtAj4p6ocxPx7C2hhebWTsWFwqSGxy2VL/jb2qyI5CRGS0rG3s8NSYt/FEiDdKdQpeX30Qv61dKDsWVSDzinxxU1VQcAM5l88CcOVS/0REFczcTItP+zVCFWstwg6OQ4e4Q4jIPI2wYZ+I7QPIuPC/aAU4H7cPv1uOww7rKahRxUZ2HCIio6fVavDmo/VhX7uFeByevAgH5o3i5oxGiMWlAmQmRIpjno0nNBqN7DhERCZBHV0JH/YxoupNEY9D09fg0KwBKC4qlB2NyhGLSwUwv1y2IFK+e7DsKEREJid0wGuIDvkYJYoWzbK349iM3ijIz5Udi8oJi0sFqJYbL472tUNlRyEiMknNHhuN+LbzUKBYIPhGBI7OehLZBcWyY1E5YHEpZ5kZV+CllK0l4Nugtew4REQmK7jjUzjb7XukKs6Ynt0VAxdF4mouTxsZOhaXcpYUV7YjdLLGE46u7rLjEBGZtKDw7sgYHoUk2/qIv5SNJxdE4NL1PNmx6CGwuJSzvHNR4pjqUF92FCIiAlDf1wOrR4fD08kaNhnxKJgdhqTTsbJjkTEUlz179qBXr17w9Cy7GmfDhg0wNLsK62FBSU9k1+wmOwoREf2pdlV7rBkdjmm2y1BHSYLDsl5IPFI2Qk6GRa+KS15eHho3bow5c+bAECmKgjUZ3pheMgjOzZ6QHYeIiG7j6WwL7+dXI8HMD87IRrV1/XBs/2bZsciQV87t3r27uN2rwsJCcbspOzsbMl28fgPX8opgYaZBkKej1CxERPR3Lu41YDFuO47N6Y36RUdRZ9uziM3/EsGdnpYdjQxxxOV+TZ8+HU5OTrdu3t7eUvMkHItGS208mnmYwcrcTGoWIiK6OwcnF9SZuBWHbVvCWlOMBntfQPRP82THIlMoLlOmTEFWVtatW3JystQ8tvE/YpnlNEzQrpCag4iI/ntzxoYTN+KgU1eYa3TIObgcS/9Q95gjfadXp4rul5WVlbjpC6drR8VR69VMdhQiIvoP5haWaDpuObYseQ+TEhvhxqYTyCwowfiO/tyuRY8Z9IiLPlH3wqhZlCjuewS1kh2HiIjugdbMDN1GvIPRnRqJxzN3nMaq7+Zyc0Y9xuJSTpJPx8JGU4QcxQZedRrKjkNERPdIHV0Z38kf7/QKwnizdRhw7nXEzH6KmzPqKb0qLrm5uYiNjRU31blz58T9pKQk6LuMhAPimGTlJxo8EREZlqGtaqFNaHOxOWPzrF8RP7MPCm5wlV19o1fFJTo6Gk2aNBE31aRJk8T9qVOnQt8pKWVlK8eZK+YSERmqZo+NQXybuWJzxib5+3FmRjfkZF2THYv0dXJu+/btxSJuhsgx84Q4mtcIlh2FiIgegrqmyzFbJ/huGy7WekmY3RnFz/8k1oAh+fRqxMVQleoUvFn4DF4vHgHX+h1kxyEioodUv2UPXO67BtfgCP/SRGTN74aUq1myYxGLS/k4fzUPMUW+WKftDJ9adWXHISKicuAf3Bo5T/+MK3DD3MJu6LcwGmfTc2XHMnksLuUg/lJZC69X3RFmWl77T0RkLHwDgqF7MQoxzj2RklWAJ+dHIP5ipuxYJo3FpRyUxG/EQLOdaO3K2edERMbGs6obVo0OR4MajkBeOpRFHXA8YovsWCaLxaUcBCSvxDSLb9DG/LjsKEREVAHc7K2wfFQYprn8goaaM6i9dTCO/MbtXWRgcXlIik4H78IEcd/Fr7nsOEREVEEcrC3QbuwCxNqEi80Zg35/AdE/L5Ady+SwuDyky0kJcEIeihQz+AQ2lR2HiIgqkLWtPepP3Ihox86w0JQiJHoyolZ+JDuWSWFxeUippyLFMcm8JiytrGXHISKiCmZhaYWQ8SsRVbUftBoFoSemI2LJZDECTxWPxeUhFSQdFsdrjvVkRyEiokqibu3SYswiRHiPEo+rn1uPj3+KgU5nmIuoGhK9WjnXENlejRdHpVrZzqJERGQaNFotwkd8hj9WVsNrsW64GJmG9KKj+PiJhjA347hAReF39iGo2xPUKCibmOtUp5nsOEREJEHrAa9gYr9OYh2vtYcuYtbXX3NzxgrE4vIQ0nIK0bHgEwwqegO+9UJlxyEiIkmeaOqF+c80RXeLQxifMhlnZnRHbvZ12bGMEovLQ66YmwV7pFcNhY2dvew4REQkUecgD4ztFoICWKF+0RFcnt0Z19Mvy45ldFhcHsKxlGxxbODpJDsKERHpgfqteuJyn1W4Dgf4lyQge15npF48IzuWUWFxeQi14mbjNfNlCHe6KjsKERHpCf8mbZH91E9IhSt8dcnQfd0VSQlHZccyGiwuD6F55haMNt+EIMdC2VGIiEiP+AaGQDdsC5I1nqiOdNj9+ChOJZySHcsosLg8IPW8ZTWki/veQZyYS0REd6ruGwDb0dtxxqw2filpgX4/nMOBc9dkxzJ4LC4P6OKJshVzL2qqw7GKq+w4RESkh1w9vFB13E78UmMicgpLMfibKOw6kSo7lkFjcXlAuecPiWOqXYDsKEREpMccnVywdEQYOgS6Q1dShJLlAxG9aaHsWAaLxeUBWaaVTbQqcm8oOwoREek5G0szLBjcFB/6HkZnbTRCDv4PUas+kR3LILG4PCD3vNPiaO8bIjsKEREZAAszLfo99xai3B4v25zx+IeIWPoaN2e8TywuDyAnNwd2urI1XGpwxVwiIrqfzRlf+AaRXiPE4/Dz8xC14AWWl/vA4vIAjqcVIaRwAfpYLoSLew3ZcYiIyMA2Zwwb+QUi674iHoelLsfB2YNQUlwkO5pBYHF5APFixVwNqnrVkR2FiIgMVNjAt3Aw+EOUKFoEXf8N7323CQXFpbJj6T1z2QEM0bFLWeLIpf6JiOhhNO8zFodtq+CLPVewN8EKiUsPYuGzzWBvxbfnf8IRlwfwdMJELLD4As0cuJAQERE9nCZdnsGYoUNhZ2mG/Weu4s15PyIz44rsWHqLxeU+3cjLQUjxYXQ1i4a/l7vsOEREZARa+rlh2agwNLW5grevv47MuZ2Qdumc7Fh6icXlPiWdOAgzjYKrcELV6r6y4xARkZFo7F0Fnw9ogmKNJWrqklG6qDMuJsbLjqV3WFzuU+bZaHG8ZF1XzAwnIiIqLzUDQ1A6bKvYTkbdnNH6h544E1e2xQyV4TvvfdJcPiKOea71ZUchIiIj3ZzR+vmyzRndkImqa/viZNSvsmPpDRaX++SSfUIcrbybyI5CRERGyq2aN9zGbscJi/pwRD58Nw9CzD6WFxWLy30oKiyAT8l5cd+jbpjsOEREZMScnN1Qc8I2HLFujqNKbQz5JQ8/H0mBqWNxuQ/nki7gmFILaXCGZ826suMQEZGRs7FzQL2Jv2BN3c+Rq7PEuBWH8WPUBZgyFpf7cCTTFn2L3sOE6ss4MZeIiCqFpZUVPhrYCoNCfaAowLWfpyJy6esmu78Rl+a7D/EpZSvm1veqIjsKERGZEDOtBh/0aYD6JScw8NgG4DwQuTAToc99ZXK/SJvWV/uQTly8Ko4NanCpfyIiqlwajQYDn+yPSL9J4nHYlR9x8MtnUFpSAlPC4nKP1B+MxWkDsMVyMhpVKZQdh4iITFTYM2/jYKP3UKpo0OL6Lzgyoy8KC/JhKlhc7tHFxKNw0NyAryYNPl4+suMQEZEJa/74eBxtOQtFijlC8vbg9IweyMvJhClgcblH6acPiGOSZR2YmXNqEBERydWk6xCc6rgY+YoVGhYexsxF3yAzvwjGjsXlHpVcihXHLKd6sqMQEREJDdv2RvJjK/EBRmFRWiD6L4hAanYBjBmLyz1yuH5MHDU1gmVHISIiuiWg6SPoP+ZteDha4XRqLp6b9wsunS17zzJGLC73QFdaCp/CBHHf1a+F7DhERER3qOvhgDWjW6K+iw4f5b0D6+964Gx8FIwRi8s9uHzhpJiYq06C8gngHkVERKR/vF1ssfTZxrAy18IVmXBb0xcnD2yHsWFxuQcJV7KxobQloqzCYWFpJTsOERHRXVWt5gPXsTtxwiIIjsiD7y8DcXT3WhgTFpd7YFbVDxtrv4v9IZ/JjkJERPSvnFzc4DthG45aN4eNpgiBu0bh0ObFMBYaRVF3PjAO2dnZcHJyQlZWFhwdHWXHISIikqaosABxXz2Fpjm7oFM0iGzyEVr2GQ1Df//miAsREZERsrSyRvD4NYhy7YMrcMYrkTaYt/sMDB2LCxERkZEyMzdHixeXYF3TH5ACN3y89SSmbzkBQz7ZwuJCRERkxDRaLcY+1gpTugeKxxf2rkDUl0MMdnNGrl1PRERkAp5vVwfu2ix02zEXNteKcGhmX9QfuxJW1rYwJBxxISIiMhF924TgRPinKFLMEJK7B6dmPor83CwYEhYXIiIiExLSbShOdfxGbM7YqCAGSTO7IutaGgwFiwsREZGJadi2L5J6LUcW7BBYcgJXv+qMjJQLMAQsLkRERCYosFlHXOu/EelwRm3deWxYMh1JV/Oh71hciIiITFStoOYoenYrfjB/Ah/m9ES/+ftx8ko29BmLCxERkQmrUTsQXcbNRUA1J6TlFGLQ/L04Ebsf+orFhYiIyMS5O1pj5XPhaObtgHdLZ8N3fR/E7VkHfcTiQkRERHCytcB3Q5vAx64EtppCBOwciUNblkDfsLgQERGRYGvngICJmxBj3x6WmlI0jpyIA2tnQJ+wuBAREdEtVlY2CJ6wFlEuj8FMo6BF3DuI/H4q9AWLCxEREf19c8ax3yKi+rPicdiZWdj79St6sTkjiwsRERHddXPG8Oe/RGTtcShQLPDVmWp4fX08SnVyywuLCxEREf2jsGffx9YOm3EA9bD8QBLGrTgMncTywuJCRERE/6pPuxb46ukQWJhpEODhAK1WA1nMpX1mIiIiMhg9G1VHQDUH1KlqJzUHiwsRERHdEz93e8iml6eK5syZg5o1a8La2hqhoaE4cOCA7EhERESkB/SuuKxcuRKTJk3C22+/jUOHDqFx48bo2rUr0tLSZEcjIiIiyTSKPlyUfRt1hKV58+b46quvxGOdTgdvb2+89NJLeO211+742MLCQnG7KTs7W3xsVlYWHB0dKz07ERER3T/1/dvJyeme3r/1asSlqKgIMTEx6NSp063ntFqteBwREfG3j58+fbr4Qm/e1NJCRERExkuviktGRgZKS0vh4eFxx/Pq4ytXrvzt46dMmSLa2c1bcnJyJaYlIiKiymbQVxVZWVmJGxEREZkGvRpxcXNzg5mZGVJTU+94Xn1crVo1abmIiIhIP+hVcbG0tETTpk2xc+fOW8+pk3PVx+Hh4VKzERERkXx6d6pIvRR6yJAhaNasGVq0aIGZM2ciLy8Pw4YNkx2NiIiIJNO74jJgwACkp6dj6tSpYkJucHAwtm7d+rcJu0RERGR69G4dl8q6DpyIiIj0g8Gu40JERET0b1hciIiIyGDo3RyXh3HzrJc65ERERESG4eb79r3MXjGq4pKTkyOOXPqfiIjIMN/H1bkuJjM5V13zJSUlBQ4ODtBoNA/8Ojc3a1S3EOAk38rB73nl4ve78vF7Xrn4/Tas77laRdTS4unpKfYoNJkRF/WL9fLyKrfXU7/x/IGvXPyeVy5+vysfv+eVi99vw/me/9dIy02cnEtEREQGg8WFiIiIDAaLy12oO06//fbb3Hm6EvF7Xrn4/a58/J5XLn6/jfd7blSTc4mIiMi4ccSFiIiIDAaLCxERERkMFhciIiIyGCwuREREZDBYXG6zZ88e9OrVS6zcp668u2HDBtmRjNr06dPRvHlzsdKxu7s7+vTpg1OnTsmOZdTmzZuHRo0a3VogKjw8HFu2bJEdy2R89NFH4t+WCRMmyI5itN555x3xPb79FhgYKDuWUbt06RKeeeYZuLq6wsbGBg0bNkR0dHSFfT4Wl9vk5eWhcePGmDNnjuwoJuH333/Hiy++iMjISGzfvh3FxcXo0qWL+O9AFUNdWVp984yJiRH/sHTo0AG9e/fGsWPHZEczegcPHsSCBQtEcaSKVb9+fVy+fPnW7Y8//pAdyWhdv34drVq1goWFhfgl6Pjx4/j888/h7OxcYZ/TqJb8f1jdu3cXN6ocW7duvePx0qVLxciL+qbatm1babmMmTqieLsPP/xQjMKo5VH9x54qRm5uLgYNGoRFixbhgw8+kB3H6Jmbm6NatWqyY5iEjz/+WOxPtGTJklvP1apVq0I/J0dcSG9kZWWJo4uLi+woJqG0tBQrVqwQI1zqKSOqOOrIYs+ePdGpUyfZUUxCQkKCOOVfu3ZtURiTkpJkRzJaP/30E5o1a4Ynn3xS/OLZpEkTUdArEkdcSG929lbP+6tDjg0aNJAdx6jFxcWJolJQUAB7e3usX78eQUFBsmMZLbUcHjp0SJwqoooXGhoqRm8DAgLEaaJ3330Xbdq0QXx8vJhPR+Xr7NmzYtR20qRJeP3118XP+bhx42BpaYkhQ4agIrC4kN78Rqr+w8Jz0RVP/Qc9NjZWjHCtWbNG/OOizjdieSl/ycnJGD9+vJjDZW1tLTuOSbj9dL86n0gtMr6+vli1ahVGjBghNZux/tLZrFkzTJs2TTxWR1zUf8vnz59fYcWFp4pIurFjx2LTpk3YtWuXmDxKFUv9TcjPzw9NmzYVV3apE9JnzZolO5ZRUudrpaWlISQkRMy7UG9qSZw9e7a4r56uo4pVpUoV1K1bF4mJibKjGKXq1av/7ZeeevXqVejpOY64kDTqNlkvvfSSOFWxe/fuCp/QRf/8G1NhYaHsGEapY8eO4tTc7YYNGyYuz508eTLMzMykZTOlidFnzpzB4MGDZUcxSq1atfrbMhanT58Wo1wVhcXlLz/gt7fyc+fOiSF1dbKoj4+P1GzGenpo2bJl2Lhxozj3fOXKFfG8k5OTWAuAyt+UKVPEULr685yTkyO+/2pp3LZtm+xoRkn9uf7rnC07Ozux3gXnclWMV155RVw9p75xpqSkiN2K1YL49NNPy45mlCZOnIiWLVuKU0X9+/fHgQMHsHDhQnGrMOru0FRm165d6k7Zf7sNGTJEdjSjdLfvtXpbsmSJ7GhGa/jw4Yqvr69iaWmpVK1aVenYsaPy66+/yo5lUtq1a6eMHz9edgyjNWDAAKV69eriZ7xGjRricWJiouxYRu3nn39WGjRooFhZWSmBgYHKwoULK/TzadT/qbhaRERERFR+ODmXiIiIDAaLCxERERkMFhciIiIyGCwuREREZDBYXIiIiMhgsLgQERGRwWBxISIiIoPB4kJEREQGg8WFiPTS+fPnodFoxC04OPihX+/ma6mb7hGR4WJxISK9tmPHDuzcufOhX+fy5cuYOXNmuWQiInlYXIhIr6kbEqq3h1WtWjWxgScRGTYWFyKqcOnp6aI4qDvI3rR//35YWlre92jK0KFD0adPH/FaHh4e4tTPe++9h5KSErz66qtiN3cvLy8sWbKkAr4SIpLNXHYAIjJ+VatWxeLFi0Xh6NKlCwICAjB48GCMHTsWHTt2vO/X++2330Q52bNnD/bt24cRI0aIItS2bVtERUVh5cqVeP7559G5c2fxcURkPDjiQkSVokePHhg1ahQGDRqE0aNHw87ODtOnT3+g11JHVWbPni0K0PDhw8UxPz8fr7/+Ovz9/TFlyhQxmvPHH3+U+9dBRHJxxIWIKs1nn32GBg0aYPXq1YiJiYGVldUDvU79+vWh1f7/713qKSP1dW8yMzMT82LS0tLKJTcR6Q+OuBBRpTlz5gxSUlKg0+nE5c4PysLC4o7H6mXOd3tO/TxEZFw44kJElaKoqAjPPPMMBgwYIE7tjBw5EnFxcXB3d5cdjYgMCEdciKhSvPHGG8jKyhJzUyZPnoy6deuK+SlERPeDxYWIKtzu3bvF4m/ff/89HB0dxfwU9f7evXsxb9482fGIyIDwVBERVbj27dujuLj4judq1qwpRmDu19KlS+9ajP7qYebQEJH+YnEhIr3WsmVLsVeRuk7Lw7C3txeL1FlbW5dbNiKqfCwuRKSX1IXjEhISxP0HvWz6drGxsbculSYiw6VRFEWRHYKIiIjoXnByLhERERkMFhciIiIyGCwuREREZDBYXIiIiMhgsLgQERGRwWBxISIiIoPB4kJEREQGg8WFiIiIYCj+D9V0o/coYK+vAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# this is here so that both are shown in the plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(solution[\"x\"](t=T), solution[\"phi\"](t=T), label=\"Uniform mesh\")\n",
    "ax.set_xlabel(\"x [m]\")\n",
    "ax.set_ylabel(r\"$\\phi$\")\n",
    "ax.legend()\n",
    "\n",
    "T = 10\n",
    "model_symbolic = TwoTankModel()\n",
    "L_min = pybamm.InputParameter(\"L_min\")\n",
    "L_mid = pybamm.InputParameter(\"L_mid\")\n",
    "L_max = pybamm.InputParameter(\"L_max\")\n",
    "geometry_symbolic = {\n",
    "    \"left domain\": {x_left: {\"min\": L_min, \"max\": L_mid}},\n",
    "    \"right domain\": {x_right: {\"min\": L_mid, \"max\": L_max}},\n",
    "}\n",
    "\n",
    "\n",
    "# mesh and discretise\n",
    "submesh_types = {\n",
    "    \"left domain\": pybamm.SymbolicUniform1DSubMesh,\n",
    "    \"right domain\": pybamm.SymbolicUniform1DSubMesh,\n",
    "}\n",
    "var_pts = {x_left: 20, x_right: 20}\n",
    "\n",
    "param = pybamm.ParameterValues({\"D_l\": 2, \"D_r\": 6})\n",
    "param.process_model(model_symbolic)\n",
    "param.process_geometry(geometry_symbolic)\n",
    "\n",
    "mesh = pybamm.Mesh(geometry_symbolic, submesh_types, var_pts)\n",
    "\n",
    "spatial_methods = {\n",
    "    \"left domain\": pybamm.FiniteVolume(),\n",
    "    \"right domain\": pybamm.FiniteVolume(),\n",
    "}\n",
    "disc = pybamm.Discretisation(mesh, spatial_methods)\n",
    "disc.process_model(model_symbolic)\n",
    "\n",
    "# solve\n",
    "solver = pybamm.IDAKLUSolver()\n",
    "\n",
    "t = np.linspace(0, T, 100)\n",
    "\n",
    "\n",
    "def solve_model_with_symbolic_mesh(L_min, L_mid, L_max):\n",
    "    t_eval = [0, T]\n",
    "    solution = solver.solve(\n",
    "        model_symbolic, t_eval, inputs={\"L_min\": L_min, \"L_mid\": L_mid, \"L_max\": L_max}\n",
    "    )\n",
    "    return solution\n",
    "\n",
    "\n",
    "solution = solve_model_with_symbolic_mesh(1, 3, 6)\n",
    "ax.plot(solution[\"x\"](t=T), solution[\"phi\"](t=T), label=\"Symbolic mesh\", linestyle=\"--\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Showing the difference in time needed to change the parameters\n",
    "\n",
    "Finally, we show the difference in time needed to change the parameters. Using the symbolic mesh results in a ~10x speedup due to the omission of the model build and discretisation steps. While this speedup is relatively small for one iteration, it becomes significant when iterating over many parameter values such as in large parameter sweeps or optimization problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for uniform mesh: 1.2151330830020015 seconds\n",
      "Time for symbolic mesh: 0.10098379100236343 seconds\n",
      "Ratio: 12.032951733546648\n"
     ]
    }
   ],
   "source": [
    "import timeit\n",
    "\n",
    "time_uniform = timeit.timeit(lambda: solve_model_with_uniform_mesh(1, 3, 5), number=100)\n",
    "time_symbolic = timeit.timeit(\n",
    "    lambda: solve_model_with_symbolic_mesh(1, 3, 5), number=100\n",
    ")\n",
    "print(f\"Time for uniform mesh: {time_uniform} seconds\")\n",
    "print(f\"Time for symbolic mesh: {time_symbolic} seconds\")\n",
    "print(f\"Ratio: {time_uniform / time_symbolic}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "sila",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
